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Abstract 


Objectives:  Permafrost  and  seasonal  ground  ice  are  key  factors  that  control  the  routing  of  water 
above  and  below  the  land  surface  in  interior  Alaska.  Hence,  frozen  ground  affects  water 
resources,  ecosystem  state,  landscape  evolution,  and  soil  stability.  Despite  its  hydrologic, 
ecologic,  and  geotechnical  importance,  the  spatial  distribution  of  permafrost  and  its  relation  to 
subsurface  water  flow  remains  poorly  understood  for  many  areas  of  interior  Alaska  largely  due 
to  its  remoteness  and  inaccessibility.  This  study  was  aimed  at  improving  the  knowledge  base  by 
addressing  the  following  objectives:  (1)  developing  a  new  numerical  modeling  tool  for 
simulation  of  groundwater/permafrost  interaction,  (2)  demonstrating  and  evaluating  geophysical 
methods  for  mapping  of  permafrost  distribution  at  several  test  sites,  and  (3)  elucidating 
permafrost-controlled  surface-water/groundwater  exchanges.  The  last  objective  was 
accomplished  using  physics-based  modeling  to  evaluate  local  settings  that  were  characterized  by 
field  efforts  and  by  evaluating  larger  scale  interactions  based  on  generic  understanding  of 
permafrost  patterns  based  on  geophysical  and  other  data. 

Technical  Approach:  (1)  The  pre-existing  U.S.  Geological  Survey  (USGS)  computer  simulation 
code  SUTRA  (Saturated-Unsaturated,  variable-density  ground-water  flow  with  solute  or  energy 
TRAnsport)  was  enhanced  to  simulate  freezing  and  thawing  of  groundwater  in  saturated  media. 
An  additional  enhancement  addressed  development  of  model  capability  to  simulate  freezing  and 
thawing  of  groundwater  in  variably  saturated  media.  These  developments  provide  a  tool  that 
enable  improved  understanding  of  the  interplay  between  frozen  ground  and  groundwater  to  be 
developed.  (2)  Complementary  ground-based  geophysical  techniques  were  used  to  map,  between 
2011  and  2014,  the  details  of  fine-scale  subsurface  geology  and  permafrost  distributions  at 
several  specific  study  sites  (10s  to  100s  of  meters)  in  the  Yukon  Flats  of  east-central  Alaska.  A 
supplemental  effort  provided  airborne  geophysical  characterization  of  large-scale  permafrost 
distribution  in  interior  Alaska  (10s  to  100s  of  kilometers).  Comparison  of  results  from  a  variety 
of  methods  enabled  multi-scale  interpretation  of  subsurface  cryologic  and  geologic  features.  (3) 
Results  from  the  small-  and  large-scale  geophysical  mapping  were  then  used  as  the  basis  for  a 
series  of  hydrologic  model  evaluations  (based  on  SUTRA  and  other  methodologies)  that 
addressed  questions  related  to  permafrost/groundwater  interaction  and  predicted  system  response 
to  past  changes  in  climatic  and  hydrologic  conditions  and  also  to  potential  future  changes  in 
these  surface  conditions. 

Results:  Completed  research  has  resulted  in  both  basic  science  and  methodological  advances.  A 
groundwater  modeling  study  of  the  Yukon  Flats  Basin  of  Alaska  demonstrated  and  quantified  the 
increase  in  baseflow  to  major  rivers  expected  to  result  from  permafrost  degradation  along  a  thaw 
sequence.  The  analysis  showed  the  transition  from  continuous  (>90%)  to  discontinuous 
(50-90%)  permafrost  coverage  to  be  a  critical  one  in  terms  of  hydrologic  change.  Fully  coupled 
fluid  flow  and  heat  transport  simulations  using  the  enhanced  SUTRA  model  quantified  rates  of 
sub-lake  talik  (i.e.,  zone  of  year-round  unfrozen  ground  that  occurs  in  permafrost  areas) 
development  for  a  variety  of  conditions  representative  of  those  in  the  Yukon  Flats.  Lake  talik 
evolution  and  cross  sectional  landscape  simulations  using  SUTRA  demonstrated  the  importance 
of  advective  heat  transport  in  accelerating  permafrost  degradation  compared  with  heat 
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conduction  alone.  Simulations  with  an  imposed  surfaee  warming  trend  illustrated  that  permafrost 
thaws  more  quiekly  in  plaees  where  there  is  aetive  groundwater  flow,  partieularly  in  groundwater 
reeharge  areas.  Rates  of  permafrost  growth  and  thaw  depend  on  eausal  and  hydrologie  faetors. 
Field-based  geophysieal  investigations  and  SUTRA  modeling  results  shed  new  light  on 
permafrost-lake  dynamics,  providing  a  physically  based  explanation  for  the  development  of  new 
permafrost  around  reeeding  lakes  in  the  presenee  of  elimate  warming.  Observations  of  new 
permafrost  in  the  margins  of  reeeding  lakes  were  shown  eonsistent  with  the  insulating  effeets  of 
eeologieal  sueeession.  SUTRA  simulations  indieated,  however,  that  such  permafrost  will  persist 
for  no  more  than  about  70  years,  under  current  projections  of  climate  warming.  Field  results  also 
provided  insight  into  the  present-day  distribution  of  permafrost  within  the  Yukon  Flats,  forming 
a  valuable  baseline  to  whieh  future  studies  ean  eompare.  Methodologieal  advanees  inelude 
demonstration  of  a  multi-seale  approaeh  to  geophysieal  mapping  of  permafrost  and  an  integrated 
assessment  of  the  eapabilities  and  limitations  of  various  geophysieal  teehniques  in  the  absenee  of 
deep  subsurfaee  information  for  ground-truthing.  Software  tools  developed  under  the  project, 
including  SUTRA  with  cold  region  eapabilities  and  FEMIC  (Frequeney-domain  EleetroMagnetic 
Inversion  Code)  for  analysis  of  large  eleetromagnetic  data  sets,  represent  researeh  produets  that 
will  enable  applieation  of  the  projeet’s  approaeh  in  other  studies. 

Dissemination  of  the  researeh  results  of  this  projeet  have  been  made  widely  available  via 
numerous  presentations  at  national/intemational  eonferenees  and  publieation  of  journal  artieles, 
and  several  additional  journal  artieles  are  planned.  Colleetively,  these  results  have  wide-ranging 
transferability  to  the  eharaeterization  and  system  behavior  of  other  permafrost-impaeted 
environments. 

Benefits:  Knowledge  gained  and  produets  developed  within  this  study  provide  a  broad 
foundation  for  future  efforts  to  evaluate,  in  detail,  the  potential  eonsequenees  of  human-indueed 
and  elimate-ehange  stresses  in  key  areas  of  Alaska.  Eurthermore,  the  hydrologie  modeling  tools 
and  approaehes  and  the  geophysieal  methods  and  evaluation  approaehes  developed  here  are 
readily  transferrable  to  a  wide  range  of  issues  related  to  permafrost  eharaeterization  and 
dynamies  affeeting  eold-region  Department  of  Defense  installations  and  operations. 

A  geophysieal  workshop  offered  as  part  of  this  projeet  presented  eomplementary  approaches  for 
subsurfaee  characterization  and  monitoring  in  eold  regions  and,  as  a  result,  provided  direet 
teehnology  transfer. 

Software  for  the  enhaneed  SUTRA  model  and  manuals  deseribing  its  development  and 
implementation  will  be  publieally  available.  Efforts  in  transitioning  SUTRA  with  its  eapabilities 
for  eold  regions  to  a  broad  user  eommunity  will  eontinue  on  well  after  the  projeet  ends,  in  the 
form  of  regular  USGS  training  workshops  and  support  provided  by  its  USGS  developers. 
Additionally,  a  software  package  was  developed  for  modeling  and  inversion  of  frequeney- 
domain  eleetromagnetie  data.  This  will  be  published  to  enable  follow-on  researeh  and  wide-scale 
applieation  of  eleetromagnetie  mapping  of  permafrost  features. 

All  software  eodes  produeed  by  this  projeet  and  all  data  eolleeted  as  part  of  this  projeet  are 
arehived  and  publieally  available.  This  open  aeeess  will  promote  usage  of  the  field-derived 
information  and  of  the  software  tools.  The  data  sets  will  serve  as  an  important  baseline  for 
subsequent  studies. 
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1.  Objectives 

The  overarching  goals  of  the  study  were  to  improve  understanding  of  climate-change  mediated 
aggradation  and  degradation  of  ground-ice  and  permafrost  distribution  in  Alaska  on  the 
groundwater  systems  that  control  ecosystems  and  their  biodiversity  and  to  develop  integrated 
approaches  for  anticipating  future  consequences. 

Our  specific  technical  objectives  addressed  within  this  report: 

1.  Develop  quantitative  software  for  numerical  modeling  of  groundwater  flow  in  aquifers 
impacted  by  freeze/thaw  processes  to  address  the  Department  of  Defense's  (DoD’s)  need  for 
enhanced  predictive  hydrologic  modeling  in  regions  impacted  by  permafrost  and  seasonal  ice. 

2.  Develop  and  implement  a  multi-method  approach  for  determining  the  permafrost 
distribution  in  specific  areas  of  strong  hydrologic  interest  via  ground-based  geophysics  and 
auxiliary  data.  The  approach  is  to  be  broadly  applicable  to  discontinuous  permafrost  regions. 

3.  Create  numerical  simulation  models  guided  and  informed  by  the  geophysical 
characterization  results,  in  order  to  gain  general  insight  into  the  interplay  of  groundwater  flow, 
ground  ice  (permafrost  and  seasonal  ice),  and  climate. 

We  hypothesized  that  permafrost  distribution  plays  an  important  role  in  the  patterns  and 
magnitude  of  groundwater  flow  and  the  distribution  of  surface  water  bodies,  and  in  turn  that 
groundwater  flow  and  surface  water  bodies  influence  permafrost  patterns  and  evolution.  An 
integrated  approach  of  permafrost  characterization  and  coupled  flow  and  energy  transport 
modeling  was  used  to  qualitatively  and  quantitatively  examine  this  hypothesis. 

We  focused  the  field  initiative  in  the  Yukon  Flats  in  interior  Alaska  due  to  (a)  the  outstanding 
potential  to  leverage  resources  and  efforts  from  a  large,  multidisciplinary  study  supported  by  the 
U.S.  Geological  Survey  (USGS),  and  (b)  the  relevance  of  this  area  to  DoD  infrastructure 
underlain  by  similar  discontinuous  permafrost  in  interior  Alaska. 
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2.  Background 


2.1  General 

Alaskan  ecosystems  are  responding  to  climate  warming  that  has  persisted  for  several  decades 
(Hinzman  et  ah,  2005).  Variations  in  wetland  and  lake  areas,  soil  moisture,  and  stream  discharge 
magnitude  and  seasonal  patterns  are  hydrologic  changes  that  directly  influence  vegetation, 
biogeochemical  cycling,  and  habitat  sustainability.  Evidence  of  lake  area  changes  in  parts  of 
Alaska  has  been  documented  by  remote  sensing  (Riordan  et  ah,  2006;  Rover  et  ah,  2012). 
Walvoord  and  Striegl  (2007)  report  historical  increases  in  baseflow  throughout  the  Yukon  River 
Basin.  The  above  mentioned  studies  propose  a  relation  between  observed  hydrologic  change  and 
permafrost  thaw.  Yet  while  other  research  has  demonstrated  that  permafrost  is  warming  and 
thawing  in  some  areas  of  Alaska  (Lachenbruch  and  Marshall,  1986;  Osterkamp  and 
Romanovsky,  1999;  Osterkamp,  2005),  the  link  between  permafrost  degradation  and  specific 
hydrologic  changes  remains  largely  speculative.  The  challenges  in  substantiating  inferred 
linkages  and  quantitatively  evaluating  permafrost  hydrology  derive  from:  (1)  the  complex  and 
coupled  interplay  of  ground  ice  evolution  and  water  flow,  (2)  general  lack  of  detailed  subsurface 
information  about  permafrost  distribution  at  multiple  scales,  and  (3)  a  non-stationary  climate, 
thus  yielding  a  non-stationary  hydrogeologic  framework. 


2.2  Hydrologic  Models  for  Simulating  Frozen  Ground  Dynamics 

Quantitative  predictive  tools  are  needed  for  understanding  processes  occurring  in  the  present-day 
subsurface  system,  for  predicting  possible  future  changes,  and  for  evaluating  management 
options  for  mitigating  undesired  ecological  changes.  However,  subsurface  ice  distribution 
depends  in  a  complex  manner  on  the  subsurface  geologic  fabric,  the  ground  surface  energy 
budget,  and  the  groundwater  flow  itself.  Numerical  modeling  codes,  available  at  the  start  of  this 
project,  did  not  function  at  the  scale  required  for  practical  hydrologic  study  and  prediction,  thus, 
new  hydrologic  computer  codes  were  needed  to  forecast  the  hydrologic  consequences  of 
permafrost  degradation  and  the  resulting  impact  to  military  and  federal  infrastructure  and 
operations.  To  provide  a  tool  with  which  to  develop  understanding  of  the  groundwater-ground 
ice  system  and  potential  feedbacks,  the  U.S.  Geological  Survey’s  SUTRA  (Saturated- 
Unsaturated,  variable-density  ground-water  flow  with  solute  or  energy  TRAnsport)  simulator 
(Voss  and  Provost,  2002)  coupling  groundwater  flow  and  heat  transport,  was  modified  to  include 
the  water  freeze/thaw  process.  To  meet  the  additional  needs  of  cold  region  model  applications, 
generalized  boundary  conditions  were  added  and  a  lake  growth/decline  package  was  coupled  to 
the  modified  version  of  SUTRA. 
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2,3  Permafrost  and  Lithologic  Characterization  Using  Geophysics 

At  the  onset  of  this  projeet,  the  aetual  permafrost  distribution  in  Alaska  as  a  whole,  ineluding 
location  of  its  top  and  bottom,  was  largely  unknown,  except  for  discrete  locations  of  limited 
extent.  A  common  but  generally  untested  assumption  (by  field  measurements)  in  Alaska  is  that 
permafrost  does  not  exist  below  major  water  bodies  such  as  rivers  and  large  lakes  (e.g.  see 
Lachenbruch  et  ah,  1962;  Williams,  1970;  Burn,  2002,  2005).  The  correctness  of  this  assumption 
is  vital  to  our  understanding  of  the  interplay  between  surface  water  bodies  and  groundwater,  and 
its  viability  underpins  any  meaningful  prediction  of  future  changes.  The  coupled 
cryologic/hydrologic  response  of  systems  to  changes  in  ground  ice  distribution/groundwater  flow 
is  key  to  assessing  the  impact  of  environmental  change  in  Alaskan  ecosystems.  Geophysical 
methods  provide  valuable  insight  into  the  spatial  and  vertical  distribution  of  permafrost  and 
relation  to  surface  water  bodies. 

Numerous  geophysical  techniques  (ground-penetrating  radar,  electromagnetics,  electrical 
resistivity,  and  seismic)  have  been  used  to  characterize  glaciers  (see  review  by  French  et  ah, 
2006)  and  the  depth  and  extent  of  permafrost  (e.g.,  Hoekstra  and  McNeill,  1973;  Daniels  et  ah, 
1976).  Compared  to  conventional  measurements  of  hydraulic  properties  (e.g.,  slug  tests, 
pumping  tests,  tracer  experiments,  and  coring)  which  require  boreholes,  geophysical  surveys 
provide  several  advantages:  (1)  expanded  spatial  coverage,  (2)  higher  resolution,  and  (3)  lower 
cost. 

Permafrost  distribution  is  a  critical  gap  in  our  general  knowledge  of  permafrost  environments 
that  limits  the  forecasting  potential  of  hydrologic  models.  For  example,  the  impact  of  permafrost 
thaw  on  DoD  lands  in  Alaska  cannot  be  adequately  addressed  without  improved  characterization 
of  the  current  distribution  of  permafrost  on  these  lands  and  cost-effective  means  of  continued 
monitoring  at  critical  locations.  To  meet  these  needs,  a  multi-method  geophysical  approach  that 
has  wide  applicability  in  characterizing  permafrost  distribution  in  the  subsurface  was  developed 
in  this  project. 


2,4  Integrated  Hydrologic  Analyses  and  Model  Applications 

The  integration  of  analytical  and  numerical  hydrologic  modeling  approaches  with  subsurface 
lithologic  and  cryologic  characterization  enables  better  understanding  of  the  influence  of 
permafrost  on  groundwater  flow  and  vice  versa.  Integrated  modeling  enables  testing  of  possible 
explanations  of  increased  baseflow  and  observed  changes  in  lake  areas  in  interior  Alaska. 
Permafrost  information  obtained  via  geophysical  methods  helps  constrain  conceptual  models  for 
both  analytical  and  numerical  approaches,  guides  numerical  model  scenario  testing,  and  provides 
a  means  of  numerical  model  calibration.  Simulations  that  incorporate  surface  temperature 
variability  are  utilized  to  address  the  effects  of  climate  on  permafrost  distribution  and  hydrologic 
processes.  To  gain  insight  into  the  interplay  of  subsurface  water  flow,  subsurface  ice,  and  surface 


5 


water  bodies,  several  modeling  studies  were  condueted  that  utilized  information  derived  from  the 
permafrost  geophysieal  characterization  results. 
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3.  Hydrologic  Model  Development 

3.1.  Materials  and  Methods 

The  USGS-SUTRA  code  (Voss  and  Provost,  2002)  that  simulates  groundwater  flow  of  variable 
density  with  either  solute  or  energy  transport,  was  originally  released  in  1984  via  a  project 
funded  by  the  U.S.  Air  Force  Engineering  and  Services  Center,  Tyndall,  A.F.B.,  Florida,  and  has 
been  continually  updated  since  then.  The  latest  public  release  of  the  SUTRA  code  was  used  a 
basis  for  incorporating  additional  capabilities,  described  below,  that  are  necessary  for  simulating 
important  subsurface  processes  in  cold  regions. 

3.1.1  Freeze/Thaw  Processes 

SUTRA  was  extended  to  incorporate  the  phase  change  between  ice  and  liquid  water.  The 
enhanced  code  simulates  dynamic  ice  formation  when  subsurface  temperatures  fall  below  a 
specified  maximum  freezing  temperature.  This  new  capability  takes  into  account  the  latent  heat 
of  freezing/thawing,  changes  in  relative  permeability  and  effective  porosity,  and  saturated- 
unsaturated  freezing.  Volumetric  changes  in  the  porous  medium  sometimes  associated  with 
freezing,  such  as  frost  heave,  are  not  represented. 

Ice  and  liquid  saturations  are  defined  as  the  fraction  of  ice  and  liquid  water  volumes, 
respectively,  in  the  total  pore  volume.  Ice  saturation  varies  over  a  specified  temperature  range 
from  0  (thus,  liquid  saturation  =  1)  at  the  maximum  freezing  temperature,  to  1  minus  a  specified 
residual  liquid  saturation  at  the  specified  minimum  freezing  temperature.  The  liquid  water 
saturation  and  the  ice  saturation  sum  to  equal  the  total  water  saturation. 

The  development  of  the  groundwater-flow  and  energy-transport  equations  is  similar  to  that  in 
previous  versions  of  SUTRA  but  is  extended  to  include  phase  change  between  liquid  water  and 
ice.  Derivation  of  the  groundwater-flow  equation  begins  with  a  transient  total-mass  balance  for 
pore  water,  both  frozen  and  unfrozen,  and  the  total-water  saturation  is  partitioned  into  liquid- 
water  and  ice  saturations.  Expansion  in  terms  of  derivatives  with  respect  to  pressure  and 
temperature  and  incorporation  of  Darcy’s  Faw  then  lead  to  a  groundwater-flow  equation  that 
accounts  for  changes  in  saturated  and  unsaturated  storage  of  liquid  water  and  ice,  flow  of  liquid 
water,  external  sources  of  water,  and  apparent  sources  of  water  associated  with  volumetric 
expansion  or  contraction  of  liquid  water  and  ice  due  to  temperature  changes.  As  mentioned 
earlier,  volumetric  changes  associated  with  phase  change  are  not  included.  Derivation  of  the 
energy-transport  equation  begins  with  a  total  thermal-energy  balance  that  includes  liquid  water, 
ice,  and  the  solid  matrix.  These  three  constituents  are  assumed  to  be  in  thermal  equilibrium  with 
each  other  at  any  given  point  in  the  model,  the  thermal  energies  of  the  constituents  are  expressed 
as  linear  functions  of  their  common  temperature.  The  latent  heat  associated  with  freeze-thaw  is 
introduced  as  the  difference  between  the  thermal  energies  of  liquid  water  and  ice  at  the  freezing 
point,  and  the  conductive-dispersive  heat  flux  is  expressed  using  an  effective,  or  “bulk,”  thermal 
conductivity  and  is  proportional  to  the  temperature  gradient.  The  first  part  of  the  derivation 
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yields  the  “divergent”  form  of  the  energy  balanee.  To  obtain  the  “eonvective”  form  solved  by 
SUTRA,  the  term  involving  the  divergence  of  the  advective  heat  flux  is  expanded,  and  the  total- 
water  mass  balance  and  the  assumption  of  an  immobile  solid  matrix  are  used  to  simplify  the 
equation.  The  result  is  a  temperature-based  thermal-energy  transport  equation  that  takes  into 
account  changes  in  storage  of  thermal  energy  in  the  liquid  water,  ice  and  solid  matrix  due  to 
changes  in  temperature  and  the  amounts  of  liquid  water  and  ice;  advection  and  diffusion- 
dispersion  of  heat;  and  external  sources  and  internal  production  of  heat. 

The  modified  SUTRA  code  allows  for  definition  of  non-linear  freezing  and  relative  permeability 
relations.  Often,  due  to  the  lack  of  established  constitutive  relationships,  liquid  water  can,  for 
simplicity,  be  assumed  to  freeze  according  to  a  linear  function  of  temperature  with  liquid 
saturation  decreasing  across  a  defined  range,  say  from  0°C  to  -1°C.  Further,  an  increase  in  ice 
saturation  can  be  assumed  to  reduce  the  effective  permeability  of  the  geologic  medium  as  a 
linear  function  between  ice  saturation  and  permeability  or  the  logic  transform  of  permeability. 
For  example,  when  assuming  the  logic,  the  relative  permeability  can  be  specified  to  have 
decreased  (log-linearly)  by  three  orders  of  magnitude  (1x10'  )  when  ice  saturation  reaches  0.80, 
a  functionality  patterned  after  nuclear  magnetic  resonance  measurements  and  ice-permeability 
analyses  by  Kleinberg  and  Griffin  (2005).  For  higher  ice  saturation,  the  logic  of  relative 
permeability  can  be  decreased  linearly  along  a  secondary  slope  toward  a  fully-restricting 
condition  (1x10'  )  when  ice  saturation  reaches  0.99.  This  coupling  of  water  movement  with 
energy  transfer  enables  the  simulation  of  both  conductive  and  advective  heat  transport  processes 
in  cold  region  environments. 

The  SUTRA  model  for  saturated-unsaturated  freezing  is  based  on  the  following  assumptions; 

1)  The  effects  of  water  vapor  are  negligible.  Only  ice  and  liquid  water  need  be  considered. 

2)  At  a  given  temperature,  the  porous  medium  contains  more  liquid  water  when  it  is  fully 
saturated  than  it  does  when  it  is  partially  saturated.  In  other  words,  partially  desaturating 
the  medium  cannot  increase  the  liquid  water  saturation  if  the  temperature  remains 
constant. 

3)  At  a  given  temperature,  the  porous  medium  contains  more  ice  when  it  is  fully  saturated 
than  it  does  when  it  is  partially  saturated.  In  other  words,  partially  desaturating  the 
medium  cannot  increase  the  ice  saturation  if  the  temperature  remains  constant. 

4)  Some  minimum  amount  of  liquid  water  always  remains,  regardless  of  the  temperature. 

5)  Given  the  preceding  assumptions,  as  little  water  as  possible  freezes. 

6)  Total-water  saturation  is  a  function  only  of  pressure  (not  temperature). 

7)  Ice  saturation  under  fully  saturated  conditions  is  a  function  only  of  temperature  (not 
pressure). 

8)  Relative  permeability  is  a  function  only  of  liquid-water  saturation. 

9)  The  formation  of  ice  does  not  deform  the  porous  medium  (i.e.,  no  volume  expansion  of 
the  matrix),  nor  is  water  attracted  to  the  freezing  front  by  decreased  pressures  (i.e.,  no 
cryosuction).  Similarly,  the  transition  from  ice  to  water  does  not  deform  the  porous 
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medium  (i.e.,  no  volume  reduction  of  the  matrix).  These  assumptions  limit  the  model  to 
hydrological  calculations  in  the  absence  of  geomechanical  feedbacks. 

The  SUTRA  saturated-unsaturated  freezing  model  was  derived  from  a  set  of  assumptions 
regarding  the  behavior  of  functions  that  determine  the  total-water,  liquid-water,  and  ice 
saturations  and  the  relative  permeability.  The  development  took  a  macroscopic  approach  in  that 
it  made  no  reference  to  physical  processes  at  the  pore  scale.  It  is  possible,  however,  to  arrive  at 
the  same  model  by  proceeding  from  a  set  of  assumptions  regarding  pore-scale  processes: 

1)  The  porous  medium  is  treated  as  a  continuum,  and  the  liquid- water  and  ice  saturations  at 
any  point  are  aggregates  of  the  states  of  individual  pores  within  a  representative 
elementary  volume  (REV)  at  that  point. 

2)  Each  pore  can  be  characterized  as  having  a  single  effective  radius  that  influences  both  its 
tendency  to  saturate  with  water  and  its  tendency  to  freeze. 

3)  An  individual  pore  is  either  completely  saturated  (fdled)  or  completely  desaturated 
(empty),  depending  only  on  pressure.  The  smaller  the  radius  of  a  pore,  the  lower  (more 
negative)  the  pressure  at  which  it  desaturates.  So,  at  a  given  pressure,  all  pores  above  a 
certain  radius  are  desaturated  (empty),  and  all  pores  below  that  radius  are  saturated 
(fdled). 

4)  An  individual  pore  is  either  completely  frozen  or  completely  unfrozen,  depending  only 
on  temperature,  the  pore  radius,  and  whether  the  pore  is  filled  or  empty.  The  smaller  the 
radius  of  a  pore,  the  lower  the  temperature  at  which  it  freezes.  So,  at  a  given 
temperature,  all  filled  pores  above  a  certain  radius  are  frozen,  and  all  pores  below  that 
radius  are  unfrozen. 

5)  The  configuration  of  liquid  water  in  the  pore  space  in  the  presence  of  ice  is  the  same  as  it 
would  be  under  ice-free  conditions  at  the  same  liquid-water  saturation. 

In  the  modified  SUTRA  model,  a  number  of  unsaturated-flow  and  ice-saturation  functions  are 
preprogrammed,  and  their  parameters  can  be  specified  through  a  SUTRA  input  dataset.  The 
preprogrammed  water-saturation  functions  are: 

•  van  Genuchten  (1980), 

•  Brooks  and  Corey  (1964), 

•  power  law,  and 

•  piecewise-linear; 

the  preprogrammed  relative-permeability  functions  are: 

•  van  Genuchten  (1980), 

•  Brooks  and  Corey  (1964), 

•  power  law,  and 

•  piecewise-linear; 
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and,  the  preprogrammed  iee-saturation  funetions  are: 

•  exponential, 

•  power  law,  and 

•  pieeewise-linear. 

In  addition,  any  user-defined  model  that  ean  be  programmed  into  the  unsaturated  eonditions 
subroutine  is  reeognized  within  the  same  framework,  so  users  may  program  their  own  funetions 
and  specify  function  parameters.  The  different  types  of  water-saturation,  relative-permeability, 
and  ice-saturation  functions  may  be  used  in  any  combination. 

3.1.2  Generalized  Boundary  Conditions 

An  additional  accomplishment  in  model  development,  exceeding  tasks  outlined  in  the  proposal, 
was  the  incorporation  of  generalized  boundary  conditions  for  SUTRA.  Two  new  types  of 
“generalized”  boundary  conditions  were  implemented  that  facilitate  simulation  of  a  wider  range 
of  hydrologic  processes  that  interact  with  the  groundwater  model,  such  as  rivers,  drains,  and 
evapotranspiration.  For  “generalized-flow”  conditions,  inflow  or  outflow  of  fluid  mass  varies 
linearly  with  pressure,  subject  to  optional  upper  and  lower  limits  on  flow  and/or  pressure.  For 
“generalized-transport”  conditions,  inflow  or  outflow  of  solute  mass  or  energy  varies  linearly 
with  concentration  or  temperature,  respectively. 

Prior  to  the  generalized  boundary  condition  implementation,  SUTRA  supported  only  the 
following  four  types  of  boundary  conditions: 

•  sources  and  sinks  of  fluid  mass, 

•  sources  and  sinks  of  solute  mass  or  thermal  energy, 

•  specified  pressures,  and 

•  specified  concentrations  or  temperatures. 

These  four  types,  which  are  still  supported  in  the  current  version  of  SUTRA,  allow  users  to 
represent  a  variety  of  physical  processes  that  interact  with  the  groundwater  model,  such  as 
pumping  or  injection  of  water,  recharge,  pressure  due  to  surface  water  of  known  depth  and 
density,  addition  of  tracer,  and  variations  in  surface  temperature. 

To  facilitate  representation  of  a  wider  range  of  physical  process,  two  additional  types  of 
boundary  conditions  have  been  added: 

•  generalized  flow  conditions  -  sources  and  sinks  of  fluid  mass  that  depend  on  pressure, 
and 

•  generalized  transport  conditions  -  sources  and  sinks  of  solute  mass  or  energy  that  depend 
on  concentration  or  temperature,  respectively. 
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The  two  new  types  of  boundary  eondition  are  "generalized"  in  the  sense  that  they  inelude  the 
four  original  types  as  speeial  eases  while  offering  additional  eapabilities.  Like  the  four  original 
types,  the  two  new  types  ean  be  made  time-dependent  using  input  fdes. 

Like  sourees  and  sinks  of  fluid  mass  and  speeified  pressures,  generalized-flow  eonditions  require 
speeifieation  of  eoneentrations  or  temperatures  assoeiated  with  inflows  to  the  model.  In  addition, 
generalized-flow  eonditions  require  speeifieation  of  eoneentrations  or  temperatures  assoeiated 
with  outflows  from  the  model.  The  latter  requirement  enables  simulation  of  proeesses  sueh  as 
evaporation,  in  whieh  fluid  may  leave  the  model  at  a  eoneentration  or  temperature  different  from 
that  of  the  groundwater. 

3.1.3  Lake  Capability 

The  lake  eapability  simulates  the  evolution  of  lakes  as  they  exehange  fluid  mass  and  solute  mass 
or  thermal  energy  with  eaeh  other  and  with  groundwater.  This  eapability  was  fully  integrated 
with  SUTRA,  thereby  allowing  lakes  to  form,  fill,  eoalesee,  drain,  empty  and  divide.  Simulated 
lakes  interaet  with  the  groundwater  flow  model  through  boundary  eonditions.  Lakes  inerease 
boundary  pressures  and  reeeive  and  supply  water  and  solute  or  energy  that  would  otherwise  have 
entered  or  exited  the  groundwater- flow  system  direetly  through  boundary-eondition  nodes.  Lakes 
ean  also  aetivate  and  deaetivate  boundary  eonditions,  giving  the  user  eonsiderable  flexibility  in 
modeling  the  effeets  of  lakes  on  the  groundwater-flow  system. 

Given  the  topography  of  the  top  surfaee  of  the  model,  SUTRA  identifies  the  lakes  that  ean  form 
by  “ponding”  of  water  on  the  top  surfaee.  Inereases  and  deereases  in  lake  stage  ean  eause  lakes 
to  eoalesee  and  divide,  leading  to  a  hierarehy  of  “parent”  and  “ehild”  lakes  that  is  automatieally 
determined  by  SUTRA  at  the  start  of  the  simulation.  On  eaeh  time  step,  SUTRA  solves  the 
groundwater  flow  equation,  eomputes  rates  of  exehange  of  water  between  groundwater  and 
lakes,  updates  lake  stages  and  eoneentrations  or  temperatures,  eomputes  rates  of  exehange  of 
solute  mass  or  thermal  energy  between  lakes  and  groundwater,  and  solves  the  groundwater 
transport  equation.  The  ealeulations  are  performed  sueh  that  water  mass  and  solute  mass  or 
thermal  energy  are  eonserved.  To  illustrate  how  the  new  funetion  determines  the  set  of  potential 
lakes,  eonsider  the  example  shown  in  Figure  3.3.1.  For  simplieity,  the  topography  in  this 
example  varies  along  only  one  grid  direetion,  effeetively  redueing  the  dimensionality  of  the 
problem. 
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node  number 


Figure  3, 1,3,1  Example  for  nine  potential  lakes  based  on  given  surface  topography. 

To  simplify  the  explanation,  the  dimensionality  of  the  problem  has  been  reduced,  and  the 
node  numbering  includes  only  surface  nodes.  Black  dots  represent  surface  nodes.  Blue  and 
green  horizontal  lines  show  the  stages  at  which  the  various  lakes  come  into  existence  by 
splitting  of  larger  lakes  as  the  stage  decreases.  Two  different  colors  are  used  to  help 
differentiate  adjacent  (“sibling”)  lakes. 


SUTRA  identifies  potential  lakes  by  starting  with  the  highest  surfaee  node  and  working  its  way 
down  through  surface  nodes  in  order  of  decreasing  elevation,  finding  the  “sills”  at  which  lakes 
divide.  In  this  example,  the  search  begins  with  the  water  stage  set  at  1.0  (the  highest  elevation, 
which  is  found  at  nodes  1  and  31).  At  this  stage,  all  the  lake  water  is  coalesced  in  one  lake, 
which  we  call  lake  1 . 

When  the  stage  drops  below  an  elevation  of  0.9  (node  6),  lake  1  divides  into  two  lakes:  lake  2  on 
the  right  (nodes  7-31)  and  lake  3  on  the  left  (nodes  1  -  5). 

The  next  division  occurs  when  the  stage  drops  below  0.8  (node  26).  Then  lake  2  divides  into 
lake  4  (nodes  8  -  25)  and  lake  5  (nodes  27  -  30).  Similar  divisions  occur  as  the  stage  drop 
through  elevations  of  0.7  (node  12)  and  0.6  (node  20). 

In  the  end,  SUTRA  identifies  the  9  potential  lakes  shown  in  Figure  3.3.1.  The  lakes  form  a 
hierarchy  of  “parents”  and  “children”  that  is  mapped  out  as  an  inverted  tree  in  Figure  3.3.2. 

The  potential  lakes  are  conceptualized  as  distinct  objects  that  can  each  hold  a  certain  volume  of 
water.  Each  surface  node  is  assumed  to  be  overlain  by  a  "cell"  that  can  collect  water,  and  the 
capacity  of  a  lake  is  the  sum  of  the  volumes  of  all  the  cells  in  that  lake.  The  rate  at  which  the  lake 
stage  rises  as  the  volume  of  water  in  the  lake  increases  is  proportional  to  the  surface  area  of  the 
lake,  which  is  the  sum  of  the  areas  associated  with  all  the  submerged  surface  nodes  in  the  lake. 
Thus,  lake  stage  is  a  piecewise-linear  function  of  volume,  with  an  increase  in  slope  occurring 
each  time  the  lake  stage  rises  high  enough  to  submerge  another  node.  If  the  volume  of  water  in  a 
lake  is  known,  its  stage  can  be  computed  using  this  function.  Each  lake  has  its  own  stage-volume 
function,  which  is  computed  automatically. 


12 


Figure  3, 1,3,2  Hierarchy  of  nine  lakes  identified  in  the  example 
problem. 

Moving  down  the  tree  as  stage  decreases,  lake  1  divides  into  lakes  2 
and  3,  lake  3  divides  into  lakes  4  and  5,  etc.  Moving  up  the  tree  as 
stage  increases,  lakes  8  and  9  coalesce  into  lake  7,  lakes  6  and  7 
coalesce  into  lake  4,  etc.  Lake  7  is  the  "parent"  of  lakes  8  and  9; 
lakes  8  and  9  are  “sihlings,”  as  they  are  both  "children"  of  lake  7, 


The  incorporation  of  the  lake  capability  into  the  version  of  SUTRA  with  freezing  capabilities 
was  viewed  as  a  critical  model  development  aspect  for  simulating  groundwater-surface  water 
interaction  among  sets  of  cold  region  lakes  that  may  be  variably  hydrologically  connected  in 
both  time  and  space  depending  on  permafrost  configuration  and  interaction  with  climate. 


3,2  Results  and  Discussion 

The  modified  version  of  SUTRA  with  freeze/thaw  capabilities  has  been  compared  with 
analytical  solutions  (Kurylyk  et  ah,  2014a)  and  has  been  utilized  by  end  users  internal  to 
(described  in  section  5)  and  external  to  (i.e.,  Ge  et  al.,  2011;  Kurylyk  et  al.,  2014b)  the  SERDP- 
supported  project.  Version  3.0  incorporates  the  generalized  boundary  conditions  and  the  lake 
evolution  capability.  Full  documentation  of  the  added  capabilities  and  instruction  on  how  to 
apply  these  capabilities  with  SUTRA  simulations  are  provided  in  Provost  and  Voss  (2014; 
USGS-approved  and  in  final  production).  Version  4.0  incorporates  the  freeze/  thaw  capability  in 
addition  to  the  added  functions  in  Version  3.0.  Full  documentation  and  user  instructions  are 
provided  in  Voss  et  al.  (in  preparation).  The  SUTRA  3.0  and  4.0  software  will  be  published 
concurrently  with  the  reports  and  will  be  publicly  available  online. 

Results  of  simulations  using  the  modified  SUTRA  code  are  detailed  in  section  5  of  this  report.  In 
addition,  results  of  idealized  test  problems  are  discussed  below. 

3.2.1.  Freeze/Thaw:  The  Frozen  Wall 

This  example  demonstrates  the  effect  of  freezing  on  the  flow  field  in  a  groundwater  system.  The 
frozen  wall  is  an  engineered  intervention  that  creates  a  region  of  low  permeability  that  shields  a 
contaminated  portion  of  the  aquifer  from  regional  groundwater  flow  that  would  otherwise  spread 
contaminants  from  the  contaminated  site  to  a  larger  area  (i.e,.  thermosyphoning  technology  and 
applications).  A  similar  “frozen  wall  problem”  was  first  proposed  by  McKenzie  et  al.  (2007). 

The  area  considered  is  a  25  m  by  15  m  area  in  which  contaminant  exists  in  a  4  m  by  4  m  square 
region  centered  along  the  south  edge  of  the  area  (Figure  3. 2. 1.1).  Regional  groundwater  flow  in 
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this  shallow,  uniform  aquifer  occurs  uniformly  from  west  to  east.  An  L-shaped  frozen  “wall”  0.5 
m  in  width  is  emplaced  by  burying  freezing  coils  through  the  entire  thickness  of  the  shallow 
aquifer.  The  wall  is  installed  just  upstream  and  along  the  lateral  side  of  the  contaminated  area,  at 
a  distance  of  0.5  m  from  the  contaminated  zone  to  avoid  entraining  contaminants  during 
emplacement.  The  wall  extends  5  m  in  the  north-south  and  5  m  in  the  east-west  directions. 


Figure  3,2, 1,1  Model  schematic  for  the  frozen  wall  example  problem. 


The  modeled  area  is  initially  at  a  temperature  of  5  °C  and  water  of  the  same  temperature  steadily 
recharges  the  area  from  one  side.  At  the  start  of  the  simulation,  a  specified  sub  freezing 
temperature  is  set  within  the  hypothetical  wall,  which  is  held  at  this  temperature  throughout  the 
simulation.  This  may  be  considered  as  one  half-section  of  a  larger  region  containing  a  series  of 
walls  with  uniform  spacing  that  are  perpendicular  to  the  regional  groundwater  flow  direction. 
Figure  3. 2. 1.2  illustrates  the  flow  field  and  temperature  distribution  approximately  25  days  after 
the  wall  freezing  is  initiated.  Groundwater  flows  around  the  frozen  wall,  which  forms  a  very  low 
permeability  barrier.  Figure  3. 2. 1.3  shows  the  corresponding  ice-saturation  distribution. 
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Figure  3,2, 1,2  Flow  field  and  temperature  distribution  after  25  days. 

White  squares  denote  model  nodes  that  represent  the  frozen  wall.  Velocity  vectors  are 
drawn  using  a  dot  to  show  location  and  line  segment  to  show  direction  and  relative 
magnitude.  Colored  shading  represents  temperature  in  °C,  as  indicated  in  the  legend. 
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Figure  3,2, 1,3  Ice-saturation  distribution  after  approximately  25  days. 

White  squares  denote  model  nodes  that  represent  the  frozen  wall.  Velocity  vectors  are 
drawn  using  a  dot  to  show  location  and  a  line  segment  to  show  direction  and  relative 
magnitude.  Colored  shading  represents  ice  saturation,  as  indicated  in  the  legend.  The  area 
that  is  not  colored  is  ice-free, 

3.2.2.  Generalized  Boundary  Conditions:  Dam  with  Vertical  Sides  (Seepage  Face) 

This  problem  involves  two-dimensional  (2D),  constant-density  groundwater  flow  through  a  IO¬ 
meter- wide  and  10-meter-high  dam  with  vertical  sides  under  steady-state  conditions.  The  water 
level  is  held  even  with  the  top  of  the  dam  on  the  upstream  side,  where  hydrostatic  conditions 
prevail,  and  water  does  not  accumulate  on  the  downstream  side,  where  water  exits  through  a 
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seepage  faee  (Figure  3.2.2. 1).  The  seepage  faee  boundary  eondition  is  implemented  using  a 
generalized  flow  eondition. 

Figure  3. 2.2. 2  shows  the  simulated,  near-steady-state  water  table  and  flow  field  after 
approximately  82  hr  of  simulated  time.  Along  the  right  vertieal  boundary,  no  flow  erosses  the 
boundary  above  the  water  table,  and  below  the  water  table  water  diseharges  through  the  seepage 
faee.  The  slight  deerease  in  the  slope  of  the  water  table  as  it  nears  the  right  vertieal  boundary  is 
due  to  the  faet  that  at  that  boundary  the  simulated  water  table  must  eoineide  exactly  with  a  node 
at  which  />  =  0  is  specified;  the  water  table  cannot  fall  between  two  nodes  at  that  boundary.  The 
intersection  of  the  water  table  with  the  right  vertical  boundary  marks  the  top  of  the  seepage  face. 
The  semi-analytical  solution  of  Polubarinova-Kochina  (1962),  which  does  not  account  for  flow 
through  the  unsaturated  zone,  predicts  a  slightly  lower  seepage  face  height  (Figure  3.2. 2.2,  black 
circle).  The  Dupuit  formula,  which  is  exact  if  it  is  assumed  that  no  flow  occurs  through  the 
unsaturated  zone,  predicts  discharge  of  4.905x10'  m  /s  per  meter  of  dam  thickness  perpendicular 
to  the  cross  section.  This  is  6  percent  less  than  the  discharge  simulated  by  SUTRA,  which  is 
5.228  xlO'  m  /s  per  meter  of  dam  thickness. 


Figure  3,2,2, 1  Boundary  conditions  and  finite-element  mesh  for  the  dam  (seepage  face) 
example. 
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Figure  3,2,2,2  Near-steady-state  SUTRA  results  for  the  seepage  face  example. 

The  thick  line  indicates  the  numerically  simulated  position  of  the  water  table  (p  =  0),  The 
intersection  of  the  thick  line  with  the  right  vertical  boundary  indicates  the  numerically 
simulated  top  of  the  seepage  face.  The  black  circle  indicates  the  top  of  the  seepage  face 
predicted  by  the  semianalytical  solution  of  Polubarinova-Kochina  (1962),  Dots  indicate 
element  centers  from  which  flow-velocity  vectors  (thin  line  segments)  emanate.  Vector 
length  is  proportional  to  velocity  magnitude.  Vectors  are  plotted  for  every  second  element 
in  each  coordinate  direction, 

3.2.3.  Lake  Capability:  3-Lake  Example 

This  problem  involves  interaction  of  3D  groundwater  flow  with  lakes  and  transport  of  solute 
between  groundwater  and  lakes.  The  model  domain  covers  a  lO-km-by-20-km  rectangular  area 
and  is  approximately  100  m  deep,  with  a  variable  surface  topography  that  includes  three 
topographic  depressions  (Figure  3.2.3. 1).  Initially,  no  lake  water  is  on  the  surface  of  the  model. 
As  topographically  driven  groundwater  flow  recharges  at  higher  elevations  and  discharges  at 
lower  elevations,  lakes  form  in  the  three  topographic  depressions.  A  source  of  solute  at  the 
model  surface  enters  the  groundwater  and  discharges  into  one  of  the  lakes.  As  the  lake  stages 
rise,  the  lakes  coalesce,  and  their  waters  mix.  Eventually,  one  large  lake  forms.  Subsequent 
withdrawal  of  lake  water  causes  the  lake  stage  to  fall,  and  the  large  lake  splits  into  the  three 
original  lakes.  At  the  end  of  the  simulation,  the  groundwater-lake  flow  system  is  near  steady 
state. 

Figure  3. 2. 3. 2  shows  areas  of  the  model  surface  covered  by  lake  water  at  a  sequence  of  times 
during  the  simulation.  Figure  3.2. 3. 3  illustrates  how  lake-water  concentrations  change  with  time. 
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Figure  3,2,3, 1  Boundary  conditions  and  finite-element  mesh  for  the  3-lake  example. 

Green  cubes  denote  specified-pressure  nodes  that  act  as  solute  sources.  Red  cube  denotes  a 
fluid  sink  node  that  withdraws  lake  water,  (SOX  vertical  exaggeration) 


Figure  3,2,3,2  Areas  of  the  model  surface  covered  by  lake  water  (blue)  in  the  lake  example. 
Groundwater  flow  recharges  at  higher  elevations  and  discharges  at  lower  elevations, 
supplying  water  to  form  lakes  in  the  topographic  depressions.  A,  By  30  weeks,  lakes  3 
(right)  and  4  (left)  are  full  and  are  spilling  over  onto  lake  5  (middle),  B,  By  60  weeks,  lakes 
4  and  5  are  coalesced  into  their  parent,  lake  2,  C,  By  100  weeks,  all  lakes  are  coalesced  into 
one  large  lake  (lake  1),  and  lake-water  withdrawal  begins,  D,  By  500  weeks,  lake  stages 
have  declined  such  that  there  are  once  again  three  distinct  lakes,  and  inflows  nearly 
balance  outflows  in  each  lake,  so  the  flow  system  is  near  steady  state,  (SOX  vertical 
exaggeration) 
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Figure  3,2,3,3  Lake-water  concentration  as  a  function  of  time  in  the  lake  example. 
Concentration  is  expressed  in  arbitrary  units  in  which  the  maximum  source  concentration 
is  1,  When  two  lakes  coalesce  (lakes  4  and  5  at  47  weeks,  and  lakes  2  and  3  at  92  weeks),  the 
water  in  the  parent  lake  is  a  mixture  of  the  concentrations  in  the  two  child  lakes.  When  a 
lake  splits  into  its  children  (lake  1  at  119  weeks,  and  lake  2  at  347  weeks),  the  concentration 
in  each  child  lake  is  initially  the  concentration  from  the  parent  lake.  Over  time,  the 
concentrations  in  the  two  child  lakes  diverge. 
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4.  Geophysical  Investigations 

Geophysical  field  investigations  were  conducted  at  two  sites  in  interior  Alaska  in  late  summer 
2011,  2012,  and  2014.  Geophysical  methods  demonstrated  included  ground  penetrating  radar, 
frequency-domain  electromagnetics,  direct-current  resistivity,  and  time-domain 
electromagnetics.  Operation  and  information  acquired  from  these  methods  are  illustrated 
schematically  in  Figure  4.0.1. 
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Figure  4,0,1  Schematic  diagram  illustrating  setups  and  data  products  for  geophysical 
methods  used  in  this  work:  (a)  ground  penetrating  radar  (GPR);  (h)  multi-frequency 
electromagnetic  induction  using  a  GEM-2  instrument;  (c)  two-dimensional  electrical 
resistivity  (2D  Resistivity)  and  continuous  resistivity  profiling  (CRP);  and  (d)  time-domain 
electromagnetics  (TDEM), 
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The  geophysical  methods  used  in  our  work  were  selected  based  on  their  complementary  and 
overlapping  spatial  resolutions  and  coverage  (Figure  4.0.2).  In  addition  to  our  ground-based 
geophysical  methods,  another  USGS  group  recently  collected  airborne  electromagnetic  (AEM) 
surveys  under  related  research  in  Alaska  (Minsley  et  ah,  2012).  We  collaborated  with  the  other 
USGS  group  to  integrate  the  AEM  information  with  the  ground-based  geophysical  information 
collected  in  this  study.  As  illustrated  in  Figure  4.0.2,  the  various  methods  provide  information 
over  different  ranges  of  depth,  or  with  varying  spatial  resolution,  or  with  different  spatial 
coverage.  For  example,  airborne  electromagnetics  provides  the  best  spatial  coverage  of  methods 
considered  but  at  low  resolution  compared  to  several  other  methods,  such  as  ground-penetrating 
radar.  In  addition,  we  also  collected  passive  seismic,  temperature,  and  frost-probe  data.  Passive 
seismic  was  conducted  as  a  test  of  a  new  method,  not  previously  (to  our  knowledge)  used  in 
permafrost  settings.  Temperature  and  frost-probe  data  were  collected  for  ground  truth. 


Figure  4,0,2  Schematic  diagram  illustrating  the  tradeoff  in  coverage  and  resolution  for 
geophysical  methods  used  in  this  work. 


4,1  Study  Area 

The  Yukon  Flats  is  located  within  the  discontinuous  permafrost  zone  in  eastern  interior  Alaska 
(Jorgenson  et  ah,  2008).  Geophysical  surveys  were  conducted  at  two  sites  in  the  Yukon  Flats; 
Beaver  Meadow  Lake  near  Fort  Yukon  north  of  the  Yukon  River,  and  the  Twelvemile  Lake 
study  area,  which  is  situated  ~18  km  (~12  miles)  southwest  of  Fort  Yukon,  and  ~12  km  south  of 
the  Yukon  River  main  channel  (Figure  4.1.1).  The  Twelvemile  Lake  location  has  been  a  USGS 
research  site,  hence  our  research  there  was  synergistic  with  past  and  concurrent  research.  The  site 
is  remote  and  accessed  primarily  by  float  plane.  The  Beaver  Meadow  site  is  close  to  Fort  Yukon, 
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Alaska  and  accessed  by  roads,  and  thus  provided  a  more  eost-effective  test  site  for  method 
demonstration  of  geophysieal  approaehes.  Additional  field  demonstrations  were  eonducted 
around  Fairbanks  during  a  training  course  put  on  by  project  Pi’s  as  a  tech- transfer  and  outreaeh 
effort  in  August  2014. 


From  a  lack  of  inlet  or  outlet,  Twelvemile  Lake  was  stagnant  for  deeades  (until  flooding  in 
2013),  yet  numerous  relict  sediment-filled  channels  surround  the  lake.  The  two  most  obvious 
and  reeently  active  channels  are  loeated  in  the  northwest  and  southeast  eorners  of  the  lake.  The 
Northwest  ehannel  meanders  and  breaks  off  into  multiple  sub-channels.  The  southeast  channel  is 
more  defined  and  evidence  suggests  that  it  is  an  old  inlet  to  Twelvemile  Lake  with  flow 
originating  from  Buddy  Lake,  2-3  km  distant.  This  drainage  feature  originates  from  grasslands 
surrounding  Buddy  Lake  and  immediately  narrows  into  a  series  of  diseonneeted  water-filled 
thermokarst  features  and  hummoeks.  Jepsen  et  al.  (2012b)  show  an  abrupt  2  m  water  table 
elevation  drop  between  Buddy  Lake  and  Twelvemile  Lake  within  this  beaded  stream.  The 
thermokarst  beaded  stream  continues  for  approximately  500  m  before  it  widens  as  it  approaehes 
Twelvemile  Lake.  The  overall  elevation  range  between  Buddy  and  Twelvemile  Lake  is  <  20 
meters.  Vegetation  is  eharacterized  by  black  spruce  transitioning  into  aspen,  low  growth  shrubs, 
and  grasses  with  increasing  proximity  to  eurrent  or  recent  water  bodies. 


Figure  4,1,1  Location  map  of  Beaver  Meadow  and  Twelvemile  study  areas. 


4,2  Materials  and  Methods 

Field  geophysieal  eampaigns  were  conducted  with  SERDP  support  in  August  of  2011,  2012,  and 
2014  and  April  of  2012.  Additionally,  our  group  eollected  field  data  in  August  2010  with  USGS 
funding.  Complementary  geophysieal  methods  were  employed,  providing  information  about 
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subsurface  conditions  at  different  depths  or  with  different  spatial  coverage.  Co-located  surveys 
were  conducted  to  compare  results  from  different  methods  and  to  evaluate  depth  of  investigation 
and  spatial  resolution.  Several  transects  from  2011  were  reoccupied  in  2012  for  repeat  surveys  to 
evaluate  changes  over  time. 

4.2.1  Resistivity 

Direct-current  electrical  resistivity  measurements  were  collected  using  a  10-channel  IRIS  Syscal 
Pro  instrument  and  inverted  using  several  different  software  packages.  The  method  can  image  to 
a  depth  of  approximately  one  fifth  the  maximum  spread  between  electrodes  (e.g.,  100  m  for  a 
0.5 -km  line).  Resistivity  increases  in  saturated  porous  media  with  freezing,  hence  the  ability  of 
the  method  to  detect  permafrost.  In  total,  5  resistivity  lines  were  collected  at  Twelvemile  (~2 
km)  and  four  at  Beaver  Meadow  (-1.5  km).  Several  lines  were  collected  at  multiple  electrode 
spacings  to  provide  multi-resolution  information.  Whereas  longer  spacing  provides  superior 
depth  of  investigation  and  coverage,  shorter  spacings  provide  superior  resolution.  Representative 
results  are  shown  in  Figure  4.2.1 . 1 . 


Iog10(resistivity,  ohm-m) 

Figure  4,2, 1,1  Resistivity  line  2  from  Twelvemile  Lake  collected  August  2011,  with  images 
for  4-m  and  1-m  electrode  spacing,  providing  multi-resolution  information  about 
subsurface  conditions. 

To  understand  the  ability  of  resistivity  to  resolve  shallow  permafrost,  we  conducted  numerical 
experiments  on  hypothetical  permafrost  configurations  and  various  survey  geometries.  SUTRA 
models  were  run  to  produce  physically  based  distributions  of  soil  moisture  and  temperature, 
which  were  converted  to  cross  sections  of  resistivity  using  the  petrophysical  relation  Archie’s 
Law.  Synthetic  resistivity  data  were  computed  using  a  finite-difference  model,  random  noise 
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(2%)  was  added  for  realism,  and  data  were  inverted  using  the  same  eodes  and  proeedures  applied 
to  field  data.  This  analysis  provides  valuable  insight  into  survey  design.  For  example  (Figure 
4.2. 1.2),  a  2-meter  eleetrode  spaeing  is  eapable  of  revealing  shallow  resistivity  variations 
undeteeted  by  inversions  of  data  eolleoted  using  a  4-meter  eleetrode  spaeing. 
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Inversion  of  forward  model  response  with  2%  noise,  2  m  electrode  spacing 


j_ I_ I_ I_ I_ I_ I_ I_ L 


0  10  20  30  40  50  60  70  80  90 

Distance  (m) 
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Figure  4,2, 1,2  Synthetic  examples  illustrating  the  resolving  power  of  2D  resistivity  using 
different  survey  geometries.  The  top  plot  shows  the  true,  hypothetical  resistivity  cross 
section  based  on  SUTRA  model  output  for  soil  moisture  and  temperature.  The  second  plot 
shows  the  resulting  resistivity  inversion,  assuming  2%  noise  and  2-meter  electrode  spacing. 
The  third  plot  shows  the  same  true,  resistivity  cross  section  at  a  larger  scale.  The  fourth 
plot  shows  the  resistivity  inversion  resulting  from  data  collected  with  a  4-meter  electrode 
spacing. 
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4.2.2  Electromagnetics 

Frequency  domain  electromagnetic  (FDEM)  induction  data  were  acquired  using  a  Geophex 
GEM-2  instrument  along  ~30  linear  km  at  Twelvemile  Lake  and  -15  linear  km  at  Beaver 
Meadow.  The  GEM-2  is  sensitive  to  earth  electrical  conductivity  to  several  meters  to  several  tens 
of  meters  depth  depending  on  site  conditions  and  data  quality.  Data  are  acquired  rapidly  relative 
to  IDEM  or  resistivity  but  require  calibration.  Cross-sectional  images  of  conductivity  can  be 
obtained  through  inversion.  Representative  results  from  Twelvemile  Lake  are  provided  in  Figure 
4.2.2.1,4.2.2.2,  and4.2.2.3. 

A  new  software  package  for  frequency-domain  electromagnetic  inversion  (FEMIC)  was 
developed  to  meet  the  data  interpretation  needs  of  our  project.  Commercially  available  software 
is  not  amenable  to  inversion  of  large,  areally  distributed  electromagnetic  data,  as  produced  in  our 
work.  The  code  includes  a  graphical  user  interface  to  inversion  and  calibration  routines.  FEMIC 
supports  2D  and  3D  inversion  with  a  variety  of  regularization  criteria.  The  code  can  take 
advantage  of  multi-core  computers  and  can  run  on  desktop  computers  or  supercomputers.  The 
code  is  currently  undergoing  internal  testing  at  USGS  prior  to  submission  for  publication  in  a 
peer-review  journal.  The  code  and  executable  will  enter  the  public  domain,  enabling  follow-on 
research  and  expansion. 
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Figure  4,2,2, 1  Comparison  of  direct-current  resistivity  results  (top)  and  multi-frequency 
electromagnetic  data  (bottom)  for  Line  1,  Twelvemile  Lake, 
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Figure  4,2,2,2  FEMIC  inversion  results  of  data  from  GEM2  transects  for  1-m  depth. 
FEMIC  results  for  this  dataset  include  2D  vertical  cross  sections  along  transects,  hut 
results  are  displayed  here  for  only  the  1-m  depth. 
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Figure  4,2,2,3  Interpolated  GEM2  inversion  results  from  FEMIC  at  the  1-m  depth.  Results 
show  strong  correlation  with  vegetation,  which  also  is  an  indicator  of  permafrost  presence. 
Low  conductivity  (blue)  in  regions  covered  hy  spruce  indicate  permafrost.  Patches  of  low 
resistivity  elsewhere  indicate  ‘islands’  of  shallow  permafrost  formed  under  more  recent 
deciduous  vegetation  growth. 

Due  to  extreme  precipitation  patterns  over  the  summer  of  2014,  the  open  meadow  and 
willlow  shrub  area  in  the  North-Western  comer  of  Twelvemile  Lake  was  flooded  with  1-2  m  of 
surface  water.  This  unexpected  condition  presented  the  opportunity  to  evalutate  the  effects  of 
periodic  flooding  on  new  permafrost  formation.  EM  data  were  collected  using  a  raft  throughout 
the  area  where  new  permafrost  had  been  observed  in  2011,  2012,  and  2013.  Preliminary  data 
interpretation  indicates  that  new  permafrost  still  exists  in  the  flooded  area  (Figure  4.2. 2.4);  in  the 
future,  SUTRA  models  will  be  adjusted  to  included  flooding  events  to  predict  the  longer-term 
effects  on  new  permafrost  development. 
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Figure  4,2,2,4  Raw  GEM2  data  indicates  electrical  conductivity  of  the  subsurface  in  the 
2014  flooded  meadow  area.  Preliminary  review  of  the  EM  data  and  physical  ground-based 
measurements  indicate  that  at  this  early  stage  of  flooding,  the  new  permafrost  still  exists  in 
the  shallow  subsurface, 

4.2.3  Time-domain  electromagnetics 

Time-domain  electromagnetic  (TDEM)  data  were  acquired  using  a  Geonics  PROTEM  57 
instrument,  capable  of  soundings  to  several  hundred  meters  depth.  Data  were  acquired  at  9  sites 
distributed  around  Beaver  Meadow,  Twelvemile  Eake,  and  between  Eort  Yukon  and  Twelvemile 
Lake.  Representative  TDEM  sounding  results  for  a  site  between  Eort  Yukon  and  Twelvemile 
Lake  are  shown  in  Eigure  4.2.3. 1,  which  indicate  permafrost  between  depths  of  several  meters  to 
about  45  meters  at  this  location. 
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Figure  4,2,3, 1  Time-domain  electromagnetic  sounding  for  site  TM05  between  Fort  Yukon 
and  Twelvemile  Lake  for  minimum  layer  inversion  (red  line),  smooth-model  inversion 
(blue  line),  compared  to  results  from  airborne  electromagnetic  sounding  (gray  line), 

4.2.4  Ground  Penetrating  Radar 

In  August  2011  we  used  a  Geophysical  Survey  Systems  Inc.  (GSSI)  SIR-3000  control  unit 
coupled  with  a  GSSI  model  5103  400  MHz  bistatic  antenna  to  survey  the  near  surface  geology  as 
well  as  the  depth  and  spatial  extent  of  the  seasonal  thaw  zone  through  various  vegetation  regimes 
at  Twelvemile  Lake  near  Fort  Yukon,  Alaska.  The  antenna  was  hand  towed  at  0.5  m  s'\  oriented 
orthogonal  to  data  collection,  and  traces  lasted  between  80-140  ns  with  1024-2048  16-bit 
samples  per  trace.  Profiles  were  recorded  with  range  gain,  bandpass  filtering,  and  stacking  to 
reduce  noise  and  improve  signal  to  noise  ratios,  particularly  for  fiat  lying  reflectors.  Survey  lines 
were  repeated  and  expanded  upon  in  April  2012  at  Twelvemile  Lake  by  towing  a  GSSI  model 
3107  100  MHz  monostatic  transceiver  coupled  with  the  SIR-3000  behind  a  snowmobile  at  2-5 
km  hr'\  The  lower  frequency  radar  data  was  used  to  reach  greater  depths  for  an  improved 
assessment  of  permafrost  extent  and  local  near  surface  geology.  The  low  frequency  system  was 
impossible  to  use  in  August  due  to  the  thick  shrub  ground  vegetation  which  hampered  travel  with 
the  larger  low  frequency  system  relative  to  the  smaller  400  MHz  antenna.  The  snow  and  frozen 
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conditions  likely  improved  signal  penetration  relative  to  depths  expeeted  to  be  reaehed  in  the 
summer  due  to  deereased  attenuation  rates  from  less  free  water.  The  last  GPR  survey  in  August 
2012  at  Twelvemile  Lake  was  eondueted  using  the  400  MHz  antenna  and  similar  parameters  to 
resurvey  profiles  colleeted  in  August  2011.  Profiles  were  also  signifieantly  expanded  to  cover  a 
remnant  and  mostly  dry  (at  the  time)  drainage  system  between  Twelvemile  Lake  and  nearby 
Buddy  Lake. 

The  400  MHz  profiles  diseussed  here  were  ground-truthed  through  frost  probing,  shallow  (2-3 
m)  sediment  eores  and  dirt  pits  with  plaeed  metal  refleetors  at  known  depths  to  determine 
relative  permittivity  and  assoeiated  radio-wave  veloeities.  Shallow  (0. 5-1.0  m)  temperatures  and 
frost  probe  depths  were  also  reeorded  at  random  loeations  relative  to  these  profiles.  Ground  truth 
of  the  low  frequeney  data  was  not  possible  at  these  remote  field  sites  however  galvanie 
resistivity  and  GEM2  eleetromagnetie  profiles  were  eolleeted  along  several  of  the  GPR  profiles 
for  eomparison  to  both  the  100  and  400  MHz  GPR  results.  Approximately  '-'15  km  of  GPR 
profiles  were  eolleeted  over  the  three  field  seasons.  Interpreted  GPR  results  are  shown  for  the 
Twelvemile  Lake  area  in  Figure  4.2.4. 1. 
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Figure  4,2,4, 1  Interpreted  permafrost  distribution  based  on  ground-penetrating  radar 
collected  in  2011-2012,  in  tbe  vicinity  of  Twelvemile  Lake, 

4.2.5  Passive  Seismic 

The  horizontal-to-vertical  seismic-noise  ratio  (HVSR)  method  is  a  novel,  non- invasive  technique 
that  can  be  used  to  rapidly  estimate  the  depth  to  bedrock.  The  field  equipment  is  compact,  light¬ 
weight  and  easy  to  transport.  The  HVSR  method  uses  a  single,  broad-band  three-component 
seismometer  to  record  ambient  seismic  noise.  The  ratio  of  the  averaged  horizontal-to-vertical 
frequency  spectrum  is  used  to  determine  the  fundamental  site  resonance  frequency,  which  can  be 
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interpreted  using  regression  equations  to  estimate  sediment  thiekness  and  depth  to  bedroek  (Ibs- 
von  Seht  and  Wohlenberg,  1999;  Parolai  et  al.,  2002).  To  the  best  of  our  knowledge,  HVSR  has 
not  been  used  previously  to  investigate  depth  to  permafrost.  We  employed  the  teehnique  at 
sparse  loeations  around  Beaver  Meadow  and  Twelvemile  Lake,  ineluding  a  line  along  whieh 
resistivity  measurements  were  made  at  Beaver  Meadow  (Figure  4.2.5. 1).  Applieation  of  HVSR 
at  the  field  loeations  was  supplemental  to  our  original  geophysieal  data  eolleetion  plan.  Results 
indicate  that  HVSR  may  prove  a  cost-effective  method  for  determination  of  permafrost  depth. 
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Figure  4,2, 5,1  Representative  passive-seismic  results  overlayed  on  resistivity  cross-sectional 
image,  showing  strong  relation  between  seismic  interface  and  depth  to  permafrost  derived 
from  frost  probe. 


4,3  Results  and  Discussion 

Geophysical  results  from  our  multi-method  investigation  have  provided  information  about 
permafrost  distribution  around  Twelvemile  Lake  across  a  range  of  scales  (Figure  4.3.1), 
demonstrating  the  effectiveness  of  our  multi-scale  geophysical  toolbox  (Figure  4.0.2).  Sub-meter 
and  meter-scale  information  from  ground-penetrating  radar,  frequency  domain  electromagnetics, 
electrical  resistivity  and  passive  seismic  indicate  has  served  to  reveal  complex  lake/permafrost 
interactions,  where  new  permafrost  was  forming  in  the  margin  of  a  receding  lake.  Using 
geophysical  information  to  help  constrain  SUTRA  modeling,  we  have  explained  how  new 
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permafrost  can  form  in  seeming  contradiction  of  climate  warming  trends  (see  Section  5.6). 
Information  at  the  scale  of  meters  to  tens  of  meters  from  resistivity,  time-domain 
electromagnetics,  and  airborne  electromagnetics  has  provided  information  about  permafrost 
depth  and  large-scale  areal  distribution  of  permafrost;  these  data  provide  a  snapshot  against 
which  future  studies  can  compare  to  infer  changes  in  permafrost  distribution  over  time,  for 
assessment  of  climate-change  impacts. 

In  addition  to  providing  basic  insight  into  permafrost  distribution  and  processes,  our  SERDP 
efforts  include  methods  development  and  examples  of  methods  demonstrations.  To  the  best  of 
our  knowledge,  our  work  includes  the  first  application  of  HVSR  for  permafrost  study,  showing 
the  potential  of  this  cost-effective  method  for  use  in  remote  areas.  Development  and  publication 
of  the  FEMIC  code  will  allow  for  follow-on  studies  and  enable  more  widespread  application  of 
FDEM  and  also  AEM.  Our  well  attended  training  course,  conducted  at  University  of  Alaska, 
Fairbanks,  in  August  2014,  supported  the  transfer  and  transition  of  our  approach  to  cooperators 
from  academia,  federal  government,  and  state  agencies.  Participants  included  engineers  and  other 
staff  from  University  of  Alaska,  Alaska  Department  of  Transportation,  Alaska  Geological 
Survey,  Department  of  Defense,  Alyeska,  and  an  Alaska  legislative  staffer. 
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Figure  4,3,1  Summary  of  geophysical  methods  useful  for  characterizing  cold  region 
features  of  interest  at  variahle  scales  and  resolution. 
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5.  Hydrologic  Analyses  and  Model  Applications 

The  hydrologic  analyses  and  model  applications  completed  as  part  of  this  project  are  interrelated, 
but  can  be  partitioned  into  seven  sub  studies  described  herein.  The  first  study  (Section  5.1),  a 
groundwater  modeling  study  of  the  overall  study  region,  utilized  an  existing  standard 
groundwater  computer  code  that  simulated  permafrost  as  a  very  low  conductivity  unit  and 
required  manual  manipulation  of  changes  in  permafrost  distribution  to  assess  the  influence  of 
permafrost  on  groundwater  flow.  The  next  study  presented  (Section  5.2)  demonstrated  the  use  of 
the  modified  SUTRA  code  to  simulate  dynamic  permafrost/groundwater  flow  relations  for  a 
generalized  nested  groundwater  flow  system.  The  next  four  sub  studies  (Sections  5. 3-5. 6) 
developed  multiple  approaches  for  evaluating  cryohydrologic  processes  at  the  primary  study 
lake,  Twelvemile  Lake,  in  the  Yukon  Flats.  In  each  case,  information  from  geophysical 
investigations  was  used  to  inform  the  model/analytical  analysis.  Finally,  section  5.7  describes 
synthesis  efforts  at  the  Yukon  Flats  regional  scale,  using  geophysical  characterization,  remote 
sensing,  flow  analysis,  and  insight  gained  through  the  course  of  this  project  to  better  understand 
the  controls  on  lake  mass  balance  and  surface  water  distribution  in  discontinuous  permafrost. 

Additional  details  supplementary  to  this  report  can  be  found  in:  Walvoord  et  al.  (2012)  (section 
5.1);  McKenzie  and  Voss  (2013)  (section  5.2);  Jepsen  et  al.  (2013a)  (section  5.3);  Wellman  et  al. 
(2013)  (section  5.4);  Minsley  et  al.,  (2014)  (section  5.5);  Briggs  et  al.  (2014)  (section  5.6),  and 
Jepsen  et  al.,  (2013b)  (section  5.7). 

5.1  Regional  Groundwater  Flow  in  the  Yukon  Flats  Basin 

Groundwater  and  the  exchange  between  groundwater  and  surface  water  bodies  are  key 
components  that  influence  vegetation  community,  species  habitat,  biogeochemical  cycling,  and 
biodiversity.  Understanding  the  role  of  permafrost  in  controlling  groundwater  flowpaths  and 
fluxes  is  central  in  studies  aimed  at  assessing  potential  climate  change  impacts  on  ecosystems. 
Basin-scale  studies  in  interior  Alaska  show  evidence  of  hydrologic  changes  hypothesized  to 
result  from  permafrost  degradation  (e.g.,  Striegl  et  al.,  2005;  Walvoord  and  Striegl,  2007).  We 
conducted  a  groundwater  modeling  study  of  the  Yukon  Flats  Basin  to  assess  the  hydrologic 
control  exerted  by  permafrost  and  taliks  (unfrozen  zones  within  permafrost),  elucidate  modes  of 
regional  groundwater  flow  for  various  spatial  permafrost  patterns,  and  evaluate  potential 
hydrologic  consequences  of  permafrost  degradation  at  the  regional  scale. 

5.1.1  Materials  and  Methods 

A  regional  three-dimensional  groundwater  flow  model  of  the  Yukon  Flats  Basin  (Figure  5. 1.1.1) 
was  generated  using  the  USGS-MODFLOW  code  (Harbaugh,  2000)  to  (1)  illustrate  the  control 
that  permafrost  exerts  on  regional  groundwater  flow  systems,  and,  (2)  reveal  potential  trends  and 
behavior  in  groundwater  flow  that  might  result  from  climate  warming.  Although  MODFLOW  is 
an  isothermal  model,  it  was  selected  for  the  initial  phase  of  hydrologic  analysis  in  this  permafrost 
basin,  because  of  its  efficiency  in  simulating  regional  scale  hydrology.  Results  from  the 
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MODFLOW  simulations  with  variable,  but  static,  permafrost  distributions  provide  a  regional 
framework  for  subsequent  finer  scale  dynamic  permafrost  hydrology  simulations  using  the 
enhanced  SUTRA  model. 

The  Yukon  Flats  Basin  is  a  large  sub-basin  within  the  Yukon  River  Basin  in  Alaska.  The  basin 
includes  an  alluvial  lowland  (90  -  353  m  in  elevation  and  occupying  24,080  km  )  transected  by 
the  Yukon  River,  a  dissected  to  rolling  marginal  upland  (30  -  150  m  higher  than  the  lowlands 
and  occupying  15,460  km  ),  and  bordering  highlands  (occupying  78,800  km  ).  The  most  recent 
map  of  generalized  permafrost  distribution  and  characteristics  in  Alaska  by  Jorgenson  et  al. 
(2008),  indicates  that  Yukon  Flats  Basin  spans  the  transition  between  continuous  permafrost 
(defined  as  >90  %  coverage)  to  the  north  of  the  Yukon  River  and  discontinuous  permafrost  (50- 
90%  coverage)  to  the  south. 

The  watershed  boundary  for  Yukon  Flat  Basin,  encompassing  an  area  of  118,340  km  ,  defines 
the  model  domain.  Lateral  grid  spacing  for  all  simulations  presented  here  is  1500  m,  and  there 
are  26  vertical  layers,  ranging  from  1  -  25  m  thick.  Modeled  basin  thickness  is  400  m.  The 
hydrogeology  of  the  basin,  based  on  descriptions  in  Williams  (1962)  and  the  State  Surficial 
Geology  Map  of  Alaska  (NFS,  1999),  and  data  from  a  deep  borehole  at  Fort  Yukon  (Clark  et  ah, 
2009)  are  here  interpreted  as  belonging  to  6  lithologic  units  with  different  hydraulic  properties: 
upper  alluvium  (sand  and  gravel  mix),  lower  alluvium  (sand  to  fine  silt  mix),  loess,  bedrock, 
mountain  alluvium  and  colluvium,  and  permafrost.  Permafrost  was  represented  in  the  regional 
model  as  hydrogeologic  units  with  very  low  hydraulic  conductivity  (e.g.,  Burt  and  Williams, 
1976;  Horiguchi  and  Miller,  1983;  McCauley  et  ah,  2002)  that  extended  to  a  maximum  depth  of 
90  m  as  influenced  by  the  geothermal  gradient  and  determined  by  borehole  information  (Clark  et 
ah,  2009).  For  the  permafrost  distribution  considered  plausible  under  present  climate  conditions, 
the  model  was  calibrated  using  basefiow  estimates  for  major  rivers  with  the  Yukon  Flats  Basin 
(streamfiow  stations  denoted  in  Figure  5. 1.1.1).  Calibration  was  accomplished  automatically, 
using  the  inverse  model  (parameter  estimation)  functionality  of  MODFLOW-2000.  Using  the 
present-day  roughly-calibrated  model  as  a  reference  or  ‘base  case’,  permafrost  distribution  was 
varied  in  subsequent  simulations,  and  the  impacts  of  permafrost  distribution  on  groundwater 
flow  were  evaluated  to  assess  sensitivity,  magnitude,  and  implications  of  differences.  A  full 
sequence  of  cases  from  full  permafrost  (no  taliks,  or  unfrozen  zones,  within  the  permafrost  layer) 
to  permafrost-free  conditions  allowed  testing  the  influence  of  permafrost  on  groundwater  fluxes 
and  fiowpaths  for  a  full  range  of  possible  conditions.  Primarily  steady-state  groundwater  flow 
was  considered,  but  transient  analyses  were  also  carried  out  to  check  the  magnitude  of  variation 
resulting  from  the  simplifying  steady-state  assumption.  Groundwater  fiowpaths  were  determined 
with  post-processing  using  MODPATH  (Pollock,  1994).  Water  budgets  were  calculated  in  a 
post-processing  step  by  ZONEBUDGET  (Harbaugh,  1990).  Using  ZONEBUDGET,  total 
horizontal  and  vertical  inflow/outfiow  was  determined  for  each  of  the  following  hydrologic 
regions  within  the  Yukon  Elats  Basin: 

(1)  supra-permafrost  aquifer  (0  m  to  1-4  m  below  surface  in  regions  overlying  permafrost). 
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(2)  taliks-shallow  aquifer  (0  m  to  4  m  below  surface  in  permafrost-free  areas,  not  directly 
overlain  by  river), 

(3)  river  zone  (0  m  to  4  m  below  surface  of  river), 

(4)  taliks-deep  aquifer  (4  m  to  90  m  below  surface  in  permafrost-free  areas), 

(5)  permafrost  (1-4  m  to  90  m  below  surface;  i.e.,  regions  containing  permafrost), 

(6)  sub-talik  aquifer  (90  m  below  surface  to  aquifer  bottom  below  taliks),  and 

(7)  sub-permafrost  aquifer  (90  m  below  surface  to  aquifer  bottom  below  permafrost). 


- 1— 


Figure  5,1, 1,1  Location  map  of  the  Yukon  Flats  Basin  (after  Williams,  1962), 

Surficial  geology  represents  substrate  above  permafrost  (after  Karlstrom  et  al,,  1964), 
Streamflow  stations  with  baseflow  information  available  are  denoted.  Numbers  apply  to 
stations  described  in  Table  1  of  Walvoord  et  al,  (2012),  Inset  shows  Yukon  Flats  Basin 
location, 

5.1.2  Results  and  Discussion 

Permafrost  was  shown  to  function  as  an  effective  confining  layer  in  the  model,  resulting  in 
relatively  high  hydraulic  head  below  the  permafrost  in  low-lying  regions  of  Yukon  Flats  Basin 
(Figure  5. 1.2.1  a,  b).  Large  modeled  regions  of  upward  gradient  across  the  permafrost  layer 
coincided  with  locations  of  concentrated  lake  distributions  in  the  low-lying  Yukon  Flats, 
suggesting  a  possible  link  between  the  existence  and  persistence  of  lowland  lakes  and  upwelling 
of  deep  groundwater.  The  confining  effect  of  permafrost  also  served  to  partition  fiowpaths  into 
shallow/young  and  deep/old  components.  Model-calculated  residence  times  for  porewater  in  the 
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supra-permafrost  aquifer  (shallow)  were  on  the  order  of  lO'  years,  eompared  with  residenee 
times  of  10^  -  10"^  years  for  porewater  in  aquifers  below  permafrost. 

Overall  water  budget  results  from  the  base  ease  indieate  a  groundwater- flow  system  that  supports 
recharge  equivalent  to  55  mm/yr,  or  11%  of  the  basin-wide  PRISM-estimated  precipitation 
(PRISM  Climate  Group,  Oregon  State  University,  http://prism.oregonstate.edu,  created  2/2000). 
The  budget  results  (reported  and  schematically  illustrated  in  Figure  5. 1.2. 2)  demonstrate  the 
important  role  of  open  taliks  (unfrozen  zones)  in  the  circulation  of  groundwater  in  permafrost- 
dominated  systems,  emphasizing  the  importance  of  understanding  and  characterizing  talik  extent 
and  morphology.  Surficial  fluxes  represent  flow  across  the  water  table  for  each  shallow  zone; 
downward  fluxes  represent  recharge  (precipitation  -  ET)  and  upward  fluxes  represent 
groundwater  discharge.  For  the  base  case,  76%  of  the  total  groundwater  recharge  through  the 
ground  surface  occurs  above  open  taliks,  which  comprise  only  11%  of  the  total  area.  Flow 
through  the  supra-permafrost  aquifer  totals  43.7  x  10  m  /d  representing  25%  of  groundwater 
flow  through  the  system.  Although  the  supra-permafrost  zone  covers  89%  of  the  basin  in  area,  it 
comprises  less  than  1%  of  the  basin  by  volume.  Thus,  the  result  that  a  relatively  large  proportion 
(25%)  of  groundwater  flows  through  the  supra-permafrost  zone  signifies  relatively  rapid 
flushing. 
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Figure  5,1,2, 1  Steady-state  head  distribution  for  the  supra-permafrost  zone  (for  base  case 
PDISCl)  (a),  and  the  sub-permafrost  zone  (b). 

Areas  of  upward  gradients  across  the  permafrost  layer  are  indicated  by  the  color  scheme  in 
the  legend  (c).  The  blowup  for  the  lowlands  (within  green  outline)  shows  the  distribution  of 
lakes  >  100  m  in  diameter. 

Base  case  baseflow  is  54.0  x  10  m  /d  representing  30%  of  the  total  groundwater  flow  in  the 
system.  This  baseflow  is  comprised  of  19%  supra-permafrost  groundwater  flow,  23%  intra-taltk 
flow,  and  58%  sub-permafrost  flow.  These  results  are  composite  for  the  entire  Yukon  Flats 
Basin;  individual  sub-basins  will  have  differing  proportions  of  supra-  to  sub-permafrost 
contributions  that  are  dependent  on  basin  characteristics. 
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To  evaluate  the  importanee  of  lowland  lake  taliks  on  regional  groundwater  flow,  results  from  the 
base  ease  (PDISCl)  were  eompared  to  the  simulation  in  whieh  lowland  lakes  >  100  m  diameter 
possess  open  taliks  (PDISCl L).  Pereent  inereases  in  flow  resulting  from  lake  taliks  (PDISCl L 
eompared  to  PDISCl)  are  given  in  parentheses  in  Figure  5. 1.2. 2.  In  summary,  lowland  lake 
taliks  eomp rising  only  -1%  of  the  model  domain,  exert  a  proportionally  large  effeet  on 
groundwater  eireulation  patterns,  inereasing  some  fluxes  by  nearly  20%.  The  greatest  effeets  of 
lake  taliks  on  lateral  groundwater  flow  are  loealized  and  are  thus,  loeally  substantial. 
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Figure  5.1.2.2  Steady-state  zone  water 
budget  results  (not  drawn  to  scale)  for 
the  base  case,  PDISCl. 

Groundwater  fluxes  reported  in  units  of 
10^  m^/d.  Values  in  parentheses 

represent  the  percent  increase  in  flux  for 
PDISCIL  (with  lake  taliks)  relative  to 
PDISCl  flux. 


Selected  results  are  presented  from  permafrost  configurations  that  may  be  considered  to  reflect 
key  stages  of  a  permafrost  degradation  sequence,  including:  (1)  complete  permafrost  (PALL),  (2) 
continuous  permafrost  coverage  -  (PCONT)  (95%  permafrost),  (3)  discontinuous  I  -  cold 
(PDISCl)  (base  case,  89%  permafrost),  (4)  discontinuous  II  -  warm  (PDISC2)  (67% 
permafrost),  and  (5)  permafrost-free  conditions  (PNONE)  (Table  5. 1.2.1  and  Table  5. 1.2. 2). 
Numerous  other  intermediate  permafrost  configurations  were  simulated,  described,  and  reported 
in  Walvoord  et  ah,  (2012).  The  greatest  hydrologic  changes  to  the  system  occur  within  the 
continuous  to  discontinuous  sections  of  the  permafrost  coverage  spectrum,  and  results  presented 
here  focus  on  this  range.  Figure  5. 1.2. 3  illustrates  the  strong  control  of  permafrost  distribution  on 
groundwater  discharge  patterns.  For  the  continuous  permafrost  case  (PCONT),  groundwater 
discharge  is  focused  in  the  permafrost- free  Yukon  River  valley  and  to  upland  river  valleys  with 
the  greatest  vertical  relief  (Figure  5. 1.2. 3a).  The  widening  of  taliks  under  the  Yukon  River  and 
opening  of  taliks  beneath  all  major  tributaries  for  PDISCl  has  the  effect  of  focusing  discharge 
through  river  taliks  in  the  lowlands  (Figure  5. 1.2. 3  b).  The  increase  in  vertical  connectivity  in  the 
system,  concurrent  with  progressive  loss  of  permafrost  (PDISCl  to  PDISC2  to  PNONE,  shown 
in  the  figure)  results  in  expanded  regions  of  groundwater  discharge  in  the  lowlands  (Figure 
5. 1.2.3  b-d). 
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Model 

Short 

Name 

Permafrost 

Coverage 

(%) 

Description 

Complete  Permafrost 

PALL 

100 

No  open  taliks 

Continuous  Permafrost 

PCONT 

95 

Narrow  open  taliks  only  beneath  2 
largest  rivers-  Yukon  and  Porcupine 
Rivers 

Base  Case  Discontinuous  I 

PDISCl 

89 

Open  taliks  beneath  all  major  rivers 

Discontinuous  I  with  Lake 

Taliks 

PDISCIL 

88 

Same  as  above  +  open  taliks  beneath 
lowland  lakes 

Discontinuous  II 

PDISC2 

67 

Open  wide  (+5  km  each  side  from 
base  case)  river  taliks  +  lake  taliks 

Permafrost  Free 

PNONE 

0 

Entirely  unfrozen  material 

Table  5,1.2, 1  List  of  select  model  simulations  and  corresponding  permafrost  distribution 


Model 

Average 

Recharge'' 

(mm/yr) 

Total 

Baseflow 

(mVd) 

Groundwater 
Discharged  as 
Baseflow  (%) 

Supra- 

PermafrosC 

Contribution 

to  Baseflow 

(%) 

Median 

Baseflow 

Travel 

Times  (yr) 

Median 

Baseflow 

Path 

Length 

(km) 

PALL 

13 

7.4  X  10^ 

18 

93 

3712 

5.72 

PCONT 

28 

2.3  X  10® 

20 

47 

114 

3.56 

PDISCl 

55 

5.4  X  10® 

30 

19 

335 

3.73 

PDISC2 

146 

1.1  X  10^ 

23 

8.5 

486 

3.67 

PNONE 

458 

2.6  X  10^ 

17 

3.7 

461 

3.26 

Value  is  calculated  by  the  model,  reflects  an  average  across  the  entire  model  domain,  and  represents 
maximum  recharge  capacity  for  the  prescribed  hydrogeologic  conditions.'’  For  cases  in  which  permafrost 
is  absent,  this  component  represents  shallow  (<4  m)  aquifer  contribution. 

Table  5,1,2.2  Model  results  for  a  permafrost  tbaw  sequence  in  the  Yukon  Flats  Basin, 
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Model  results  showed  total  groundwater  flow  inereasing  with  decreasing  permafrost  coverage. 
Average  basin-wide  recharge  values  (and  thus,  total  groundwater  flow  for  steady  state)  increased 
nearly  an  order  of  magnitude  from  the  base  case  (PDISCl)  to  permafrost- free  conditions 
(PNONE)  (Table  5. 1.2.1).  Baseflow  magnitudes  also  increased  substantially  with  decreasing 
permafrost  coverage,  exceeding  an  order  of  magnitude  increase  from  complete  permafrost 
(PALL)  to  permafrost- free  conditions  (PNONE)  (Ligure  5. 1.2.4;  Table  5. 1.2.2).  Pertinent 
substantial  baseflow  change  applicable  to  the  near-term  (10  -  10  yr),  as  a  result  of  current 
warming  trends,  may  be  represented  by  the  change  between  one  of  the  stages  depicted  in  Ligure 
5. 1.2.4  and  the  following  stage.  Model  results  that  showed  an  overall  trend  in  increased  baseflow 
with  diminishing  permafrost  coverage  are  consistent  with  the  hypothesis  that  permafrost 
degradation  has  led  to  the  observed  increases  in  groundwater  input  to  rivers  in  northern  basins,  as 
noted  by  Smith  et  al.  (2007),  Walvoord  and  Striegl  (2007),  and  St.  Jacques  and  Sauchyn  (2009) 
among  others.  In  addition  to  considering  changes  in  baseflow  magnitude,  it  is  important  to 
address  the  sources  and  flowpaths  of  baseflow  and  how  these  may  change  with  permafrost 
degradation.  This  is  particularly  important  in  evaluating  potential  changes  in  aquatic  chemical 
exports  as  a  consequence  of  permafrost  thaw  (Lrey  and  McClelland,  2009).  Groundwater 
modeling  results  for  the  Yukon  Llats  Basin  indicate  a  sharp  decline  in  the  proportional 
contribution  of  lateral  flow  from  the  supra-permafrost  layer  to  baseflow  near  the  continuous  to 
discontinuous  permafrost  transition  (Ligure  5. 1.2. 5;  also  reported  in  Table  5. 1.2.1).  Changes  in 
the  magnitudes  and  ratios  of  organic  to  inorganic  aquatic  chemical  exports  would  be  expected  to 
accompany  the  changes  in  the  proportional  contributions  of  supra-permafrost  flow  to  baseflow 
illustrated  in  Ligure  5. 1.2. 5  (O’Donnell  et  ah,  2012). 
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c)  Discontinuous  II  -  warm  (PDISC2) 


d)  Permafrost  free  (PNONE) 


a)  Continuous  permafrost  (PCONT) 


1  -2  .5  b)  Base  case  (PDISCl) 


Figure  5, 1,2 ,3  Model  results  for  groundwater  discharge  patterns  with  varying  permafrost 
distributions. 
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Figure  5,1,2,4  Model-calculated  baseflow  for  selected  river  reaches  with  a  variety  of 
permafrost  distributions. 

The  streamflow  stations  are  noted  by  number  in  Figure  5,1, 1,1  map  as  follows:  Sheenjek 
(2),  Chandalar  (3),  Upper  Beaver  (6),  Lower  Beaver  (7),  Birch  (8),  and  Yukon  (12), 


■  Continuous  Permafrost  (PCONT) 

□  Base  case  (PDISCl) 

□  observed 

□  Discontinuous  II  (PDISC2) 

□  Sporadic  I  (PSPORl) 

■  Permafrost  Free  (PNONE) 
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Figure  5,1,2, 5  Model-calculated  percentage  of  supra-permafrost  (or  shallow  aquifer  in 
thawed  cases)  contribution  to  baseflow. 

Range  associated  with  thaw  patterns  beginning  in  late  discontinuous  stages  indicated  by 
thicker  line. 
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Model-Predicted  Baseflow  to  Upstream  River  Segment  (m  /day) 
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‘‘station  numbers  correspond  to  streamflow  stations  located  in  Eigure  5. 1.2.1. 
Table  5,1,2.2  Model-predicted  baseflow  for  select  permafrost  distributions. 


In  summary,  model  simulations  that  represented  an  assumed  permafrost  thaw  sequence  revealed 
the  following  trends  with  decreasing  permafrost:  1)  increased  groundwater  discharge  to  rivers, 
consistent  with  historical  trends  in  baseflow  observations  in  the  Yukon  River  Basin,  2)  potential 
for  increased  overall  groundwater  flux,  3)  increased  extent  of  groundwater  discharge  in 
lowlands,  and,  4)  decreased  proportion  of  supra-permafrost  (shallow)  groundwater  contribution 
to  total  baseflow.  These  trends  directly  affect  the  chemical  composition  and  residence  time  of 
riverine  exports,  the  state  of  groundwater-influenced  lakes  and  wetlands,  seasonal  river-ice 
thickness,  and  stream  temperatures.  Presently,  the  Yukon  Plats  Basin  is  coarsely  mapped  as 
spanning  the  continuous-discontinuous  permafrost  transition  that  model  analysis  shows  to  be  a 
critical  hydrogeologic  threshold;  thus,  the  Yukon  Plats  Basin  may  be  on  the  verge  of  major 
hydrologic  change  should  the  current  permafrost  extent  decrease.  This  possibility  underscores 
the  need  for  improved  characterization  of  permafrost  and  other  hydrogeologic  information  in  the 
study  area  via  geophysical  techniques  and  ground-based  observations. 
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5.2  Generalized  Permafrost  Thaw/  Groundwater  Flow  Interaction 

The  modified  SUTRA  eode  deseribed  in  Seetion  3.0  was  applied  to  a  nested  groundwater  flow 
system  subjeeted  to  mean  annual  sub-freezing  temperatures  as  a  means  of  evaluating  the 
importance  of  groundwater  flow  in  accelerating  climate-driven  permafrost  thaw.  This  modeling 
study  was  not  intended  to  be  comprehensive  regarding  behavior  of  ground  ice  in  groundwater 
flow  settings;  rather,  was  intended  to  identify  some  of  the  most  important  processes  and 
phenomena  that  warrant  further  study. 

5.2.1  Materials  and  Methods 

In  order  to  investigate  the  interplay  of  heat  conduction  and  heat  advection  via  groundwater  flow 
with  both  seasonal  ground  ice  and  permafrost  in  arctic  hydrologic  systems,  an  idealized  typical 
cold-regions  terrain  was  modeled  with  the  modifled  SUTRA  code.  The  region  represents  an  area 
with  undulating  topography.  For  a  similar  terrain  and  within  a  set  of  nested  groundwater  flow 
systems,  the  study  evaluated  the  impacts  of  a  key  (though  not  exhaustive)  set  of  surface  and 
hydrogeologic  controls  on  the  timing  and  pattern  of  permafrost  thaw  during  climate  warming. 
The  controls  considered  (though  only  select  results  are  discussed  in  this  report)  were: 

•  frequency  and  amplitude  of  topographic  (i.e.  water  table)  undulation, 

•  permeability  value  (‘zero’  permeability  giving  conduction-only  heat  transport,  and  also  a 
range  of  low  to  high  realistic  values), 

•  vertical  anisotropy  in  permeability  (ratio  of  horizontal  to  vertical  permeability,  ranging 
from  no  anisotropy  to  high  values  representing  the  upscaled  effect  of  low-permeability 
layers), 

•  heterogeneity  in  spatial  distribution  of  permeability  (patchiness  consisting  of  high-  and 
low-permeability  structures), 

•  lakes/streams  in  valleys  with  perennially-unfrozen  bottom  water  (none,  one,  or  more 
water  bodies), 

•  initial  permafrost  state  (permafrost  initially  much  colder  than  thaw  temperature  range, 
and,  permafrost  initially  ready  to  thaw  with  only  minor  temperature  increase), 

•  extent  of  groundwater  flow  through  permafrost  body  (no  throughflow,  and,  some 
throughflow), 

•  climate  warming  rate  (low  and  high  estimates),  and, 

•  amplitude  of  seasonal  temperature  variation  (none  and  a  range  from  low  to  high) 
superimposed  on  a  linear  temperature  trend  representing  a  warming  climate. 

The  two-dimensional,  ‘Tothian  hills’  cross-sectional  model  domain  was  created  in  the  spirit  of 
the  classic  nested  groundwater  flow  analysis  of  Toth  (1963).  The  water  table  coincides  with  the 
top  of  the  model.  The  uppermost  one  meter  of  the  domain  was  used  to  explicitly  model  heat 
transfer  through  an  equivalent  thermal  boundary  layer  that  represents  the  roughness  of  ground 
surface,  vegetation  and  snow  pack,  as  described  below.  The  mean  ground  surface  elevation  of 
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the  domain  rises  linearly  from  left  to  right,  with  a  single  uniform  slope.  Sinusoidal  hills  of 
uniform  frequeney  and  amplitude  were  superimposed  on  the  sloping  surfaee.  Figure  5.2.1 
illustrates  the  finite-element  mesh  and  the  overall  extent  of  the  domain  (5  km  length, 
approximately  2  km  depth)  for  one  example  of  topographic  undulation  superimposed  on  the 
linearly-increasing  mean  ground-surface  elevation.  The  same  value  of  the  mean  slope  (20  m/km) 
was  used  for  all  analyses.  Depending  on  the  topographic  shape  and  permeability  values,  a  nested 
combination  of  local  flow  systems  (e.g.,  one  below  each  hill)  and  a  larger-scale  domain- wide 
flow  system  (driven  by  the  regional  slope)  may  exist  in  the  domain.  Spatial  discretization 
consists  of  1  m  by  50  m  finite  elements  in  the  top  (finest)  band  of  elements,  10  m  by  50  m  in  the 
middle  band,  and  approximately  50  m  by  50  m  in  the  deepest  band.  For  some  simulations  with 
high  permeability  or  heterogeneous  permeability  a  finer  discretization  was  employed  (as  little  as 
0.5  m  vertical  spacing  in  the  upper  band  of  elements). 


Figure  5.2,1  ‘Tothian  Hills’ 
model  domain  (no  vertical 
exaggeration)  for  reference 
topographic  shape  with  finite- 
element  mesh  (coarsest  of 
several  used  for  simulations) 
and  boundary  conditions. 


Boundary  conditions  for  the  energy  balance  were  defined  as  follows.  A  no-flux  (insulated) 
condition  was  specified  along  the  vertical  sides  of  the  domain.  Along  the  bottom,  a  heat  source 
of  0.085  J/(s  m  )  was  specified,  equivalent  to  the  heat  flux  that  occurs  for  a  typical  vertical  lapse 
rate  of  26.5  °C/km  when  a  typical  thermal  conductivity  of  geologic  fabrics  of  3.21  J/(s  m  °C)  is 
employed  in  the  model.  A  specified  temperature  boundary  condition  with  a  time-varying 
temperature  value  representing  air  temperature  was  specified  all  along  the  top  of  the  domain  on 
the  upper  surface  of  the  thermal  boundary  layer. 

A  general  function  was  specified  for  air  temperature,  allowing  simulation  of  a  constant  or  linear 
increase  in  mean  yearly  air  temperature  with  time,  or  additionally  with  a  superimposed 
sinusoidal  seasonal  variation.  This  air-temperature  function,  with  increasing  mean  air 
temperature,  was  applied  to  all  simulations  in  this  study.  At  a  given  time,  t  (in  years),  the  air 
temperature,  Taiv  (in  degrees  C),  is  calculated  as: 

'^Air  ~  '^Initial  '^Slope  ^  '^Amplitude  siTl(7,Ut), 
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where  the  angle  of  the  sin  is  in  radians,  TinuM  (°C)  is  the  starting  temperature,  Tsiope  (°C)  is  the 
change  in  mean  annual  air  temperature  per  year,  and  TAmpiuude  (°C)  is  the  amplitude  of  yearly  air 
temperature  variation  (i.e.  half  of  the  difference  between  maximum  and  minimum  temperature  in 
each  one-year  cycle). 

For  all  analyses  (except  one  considering  initially  warmer  permafrost),  the  initial  atmospheric 
temperature  was  -6  °C  and  the  increase  of  mean  annual  air  temperature  was  either  +1  °C  per  100 
years  (pertains  to  results  presented  here)  or  +6  °C  per  100  years  (additional  examples  provided  in 
McKenzie  and  Voss,  2013).  These  rates  of  increase  bracket  the  very  low  end  and  very  high  end 
of  the  IPCC  (2007)  and  SNAP  (Scenarios  Network  for  Alaska  and  Arctic  Planning,  2012) 
warming  prediction  range.  The  amplitudes  of  the  yearly  temperature  variations  considered  were 
±5°C,  ±10°C,  ±15°C  and  ±20°C.  An  initial  condition  with  initial  permafrost  distribution  was 
obtained  for  a  system  dominated  by  heat  conduction  with  a  constant  sub-zero  air  temperature. 

5.2.2  Results  and  Discussion 

The  divergent  evolutions  of  permafrost  thawing  during  climate  warming  in  conduction-only  and 
advection-influenced  hydrogeologic  systems  are  described  in  reference  to  results  displayed  in 
Figure  5.2.2. 1.  In  this  simulation,  the  air  temperature  warms  from  -6  to  0  °C  during  the  first  600 
years.  Initially,  there  is  a  uniform  layer  of  continuous  permafrost  from  the  ground  surface 
downwards,  ranging  in  thickness  from  220  m  below  valleys  to  300  m  below  hilltops,  with  a 
mildly  undulating  bottom  (initial  ice  saturation  at  -600  years  in  Figure  5.2.2. 1  a). 

For  the  case  of  advection-influenced  thaw  (simulation  with  groundwater  flow),  permafrost 
completely  disappears  by  770  years  (Figure  5.2.2. 1).  This  is  about  one -third  less  than  the  time 
required  for  complete  permafrost  loss  in  the  equivalent  conduction-only  system  (1100  years).  It 
is  notable  that  the  configuration  of  residual  permafrost  is  quite  different  between  the  two  cases. 
For  conduction-only  thaw  (no  groundwater  flow),  the  residual  permafrost  bodies  are  located 
below  hilltops;  in  contrast,  for  advection-influenced  thaw,  the  residual  permafrost  bodies  are 
located  below  valley  bottoms.  For  conduction-only  thaw,  the  thinnest  parts  of  the  permafrost 
layer  thaw  through  first.  For  advection-influenced  thaw,  most  thaw  occurs  where  warm  recharge 
water  is  directed  against  the  permafrost,  directly  below  hilltops;  less  thaw  occurs  laterally 
because  groundwater  temperature  decreases  as  it  flows  towards  valley  discharge  points  due  to 
the  loss  of  latent  heat  during  thawing  that  occurred  below  the  hilltops.  Results  of  this  simulation 
analysis  demonstrate  that  where  groundwater  flows  readily  in  unfrozen  ground  surrounding 
permafrost  bodies,  groundwater  flow  is  an  important  control  on  the  temporal  and  spatial 
evolution  of  permafrost  during  periods  of  thaw.  Systems  with  groundwater  flow  will  experience 
accelerated  rates  of  permafrost  degradation  overall,  and  regional  patterns  of  residual  permafrost 
can  be  very  different  in  conduction-only  and  advection-influenced  groundwater  systems. 
Conversely,  where  groundwater  flow  in  unfrozen  ground  is  minimal,  such  as  where  the 
subsurface  has  very  low  permeability  inhibiting  flow,  where  the  subsurface  is  highly  stratified 
with  low  vertical  permeability,  or  where  ground  ice  is  spatially  continuous  over  large  areas 
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(blocking  recharge),  heat  eonduetion  is  the  primary  thaw  proeess  affeeting  permafrost  evolution 
during  elimate  warming. 
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Figure  5,2,2, 1  Comparison  of  conduction-dominated  (left)  and  advection-influenced 
(middle)  permafrost  thaw  during  climate  warming. 

Results  illustrate:  (a)  permafrost  distribution  and  groundwater  flow  field  in  the  uppermost 
0,5  km  of  the  2  km  deep  model  domain  and  (b)  temperature  distribution  in  entire  model 
domain.  Horizontal  permeability,  10'^^  m^;  vertical  anisotropy,  10, 

For  given  landforms  with  particular  hydraulic  driving  forces  for  groundwater  flow,  the 
magnitudes  of  fluid  flux  and  potential  advective  heat  transport  are  dependent  on  the  permeability 
of  the  geologic  fabric.  It  follows  then,  that  permeability  and  permeability  distribution  exert 
important  controls  on  thawing  of  ground  ice.  Additional  simulations  were  designed  and 
performed  to  quantify  the  reduction  in  permafrost  thaw  elicited  from  increases  in  permeability. 
Impacts  on  permafrost  thaw  of  two  types  of  permeability  heterogeneity  in  the  geologic  fabric, 
simple  layers  and  patchiness,  are  illustrated  in  Figure  5. 2. 2.2  and  5. 2.2. 3. 

Comparison  of  the  thaw  evolution  and  timing  from  the  three  simple  layer  heterogeneity 
situations  with  initially  continuous  permafrost  reveals  that,  if  the  high-permeability  zone  is 
located  anywhere  below  the  permafrost,  it  has  little  or  no  effect  on  the  thaw  process;  the  bottom 
and  below-permafrost  results  (Figure  5. 2. 2.2  a,  b)  are  practically  equivalent.  High  flow  in  the 
zone  (see  Figure  5. 2. 2.2  b  at  700  years)  begins  only  once  permafrost  becomes  discontinuous,  and 
by  that  time,  most  thaw  has  already  occurred,  so  the  higher  groundwater  flux  has  little  impact.  In 
contrast,  a  high-permeability  zone  located  at  the  top  of  permafrost  (Figure  5. 2.2. 2  c)  enhances 
supra-permafrost  groundwater  flow,  accelerating  thawing  from  the  top  and  causing  much  earlier 
breakthrough  and  complete  permafrost  disappearance. 
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a)  High  Permeability  Zone  Near  Bottom 


b)  High  Permeability  Zone  In  Middle 


c)  High  Permeability  Zone  At  Top* 
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Figure  5,2,2,2  Effect  of  heterogeneity  in  permeability  on  evolution  of  permafrost  thaw 
patterns:  high-permeability  layers. 


An  example  of  patehy  permeability  heterogeneity  was  also  simulated  as  shown  in  Figure  5.2. 2. 3. 
The  times  to  first  breakthrough  and  complete  disappearance  of  permafrost  in  the  patchy 
heterogeneous  case  are  200  years  and  750  years,  respectively,  and  in  the  homogeneous  case  are 
600  and  770  years,  respectively.  A  through-going  talik  forms  in  the  heterogeneous  case  in  one- 
third  of  the  time  required  in  the  equivalent  homogeneous  case,  but  complete  loss  of  permafrost 
takes  about  the  same  amount  of  time.  Inspection  of  the  patterns  of  residual  permafrost  (Figure 
5. 2. 2. 3)  indicates  that  thaw  occurs  preferentially  in  the  high-permeability  patches,  with  much 
earlier  creation  of  taliks  within  and  near  these  zones  due  to  higher  localized  groundwater  flow. 
In  contrast,  the  low-permeability  patches  preserve  permafrost  for  an  extended  time  due  to  very 
low  flow  and  essentially  conduction-only  conditions  within  the  patch.  Permafrost  is  also 
preserved  in  some  regions  below  the  low-permeability  zone  within  more-permeable  fabric, 
perhaps  as  a  result  of  the  patch  shielding  its  underside  from  groundwater  flow. 
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Figure  5,2,2,3  Effect  of  heterogeneity  in  permeability  on  evolution  of  permafrost  thaw 
pattern:  example  with  high-  and  low-permeability  patches. 


Surface  water  bodies  such  as  lakes  and  streams  have  a  thermal  impact  on  the  subsurface  region 
below  and  immediately  surrounding  the  water  body.  Most  surfaee  water  bodies,  partieularly  for 
those  that  do  not  freeze  to  the  bottom  in  winter,  support  an  iee-free  zone  (talik)  below  the  surfaee 
water  body.  Here,  the  interaetion  of  surface  water  bodies  and  nested  groundwater  flow  systems 
were  examined  by  introducing  lakes  with  variable  eonfigurations  into  the  cross-seetional 
simulations.  Simulated  thaw  evolutions  with  surface-water  bodies  in  different  valleys  are 
eompared  in  Figure  5. 2.2. 4,  and  these  may  all  be  compared  with  the  equivalent  hydrogeologie 
situation  without  surfaee  water  (seeond  panel  of  Figure  5.2.2. 1).  One  result  of  a  lake  is  the 
absenee  of  an  initial  permafrost  layer  below  its  valley.  Thus,  for  a  sub-lake  valley,  there  will  be 
no  long-lasting  residual  permafrost,  as  oeeurs  below  a  lake-free  valley  when  all  other  permafrost 
has  thawed.  Times  for  total  thaw  should  thereby  be  reduced.  This  is  borne  out  by  the  results,  in 
whieh  total  thaw  oeeurs  by  600  years  or  700  years,  depending  on  the  number  of  lakes,  in 
comparison  with  total  thaw  time  of  770  years  for  the  no-lake  situation.  Further,  it  might  be 
expeeted  that  when  there  are  two  or  more  lakes  with  through-going  taliks,  some  regional 
groundwater  flow  occurs  even  at  the  beginning  of  the  time  period,  shortening  thaw  times  due  to 
adveetion.  In  this  light,  it  may  be  noted  that  the  cases  with  a  single  lake  experienee  thaw 
(eomplete  in  about  700  years)  in  a  manner  similar  to  the  no-lake  ease,  beeause  no  regional  flow 
oeeurs  until  breakthrough  in  an  additional  location.  Note  in  Figure  5.2.2. 4  (clearest  at  early 
times)  that  flow  enters  higher  lakes,  passes  laterally  downhill  below  permafrost  bodies,  and 
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discharges  to  lower  lakes.  This  aceelerates  thaw  somewhat,  with  eomplete  thaw  oceurring  by  600 
years.  However,  despite  the  eontribution  of  large-seale  deep  groundwater  flow,  the  primary  thaw 
meehanism  in  eases  with  surface  water  remains  heat  adveetion  in  the  downward  flow  below 
hilltops. 
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Figure  5,2,2,4  Effect  of  lakes  on  evolution  of  permafrost  thaw  patterns. 


Results  of  this  modeling  study  (deseribed  more  fully  in  MeKenzie  and  Voss,  2013)  demonstrate 
how  speeified  climate  warming  rates  affeet  permafrost  thaw.  Simulated  times  for  total  thaw 
under  the  1  °C  and  6  °C  per  100  years  warming  rates  are  respectively:  1100  years  and  700  years 
for  eonduetion-only,  and  770  years  and  500  years  for  the  adveetion-influenced  ease. 


5,3,  Lake  Mass  Balance  Analysis 

Many  lakes  in  northern  high-latitudes  have  undergone  substantial  changes  in  surface  area  over 
the  last  four  decades.  In  the  discontinuous  permafrost  of  Yukon  Flats,  interior  Alaska,  these 
changes  have  been  non-uniform  across  adjacent  watersheds  (Rover  et  ah,  2012),  suggesting  local 
controls  on  lake  water  budgets.  In  the  lake  mass  balance  analysis  conducted  here,  multiple 
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mechanisms  were  explored  to  explain  the  decreasing  mass  of  the  primary  study  lake,  Twelvemile 
Lake,  in  the  Yukon  Flats  (Figure  5.3.0. 1).  Twelvemile  Lake  was  chosen  because  of  (1)  its 
substantial  reduction  (~60%)  in  surface  area  since  the  early  1980’s  at  a  rate  of  approximately  12 
cm  yh'  (p  <  0.001),  and  (2)  the  nonuniformity  displayed  between  Twelvemile  Lake  and  a 
neighboring  lake  (Buddy  Lake),  2  km  to  the  southeast,  which  has  shown  no  apparent  change  in 
surface  area  since  the  early  1950’s  (Figure  5. 3. 0.2).  Changes  in  snowmelt  mass  and  infiltration, 
permafrost  distribution,  and  climate  warming  were  considered  in  this  quantitative  analysis.  The 
analysis  was  carried  out  by  testing  the  sensitivity  of  lake  volume  to  assumed  changes  in  different 
plausible  water  fiowpaths.  The  effect  of  climate,  between  1950  and  2010,  on  the  mechanisms 
considered  are  both  direct  (evaporation,  snow/rain  partitioning)  and  indirect  (configuration  of 
permafrost). 
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Figure  5,3, 0,1  Study  site  location. 

Topographic  map  of  the  Twelvemile  Lake  watershed,  based  on  a  2,5  m-resolution  LIDAR 
DEM  (Gesch  2007)  (a).  Lake  boundary  of  Buddy  Lake  during  1984  is  similar  to  that  in 
2010  and  thus  not  shown.  Schematic  cross-section  view  of  soil  textures  and  permafrost 
distribution,  triangles  denote  lake  margin  locations  (b),  line  AA’  shown  in  panel  (a). 
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Figure  5,3,0,2  Images  of  Twelvemile  and  Buddy  Lakes  between  1952  and  2011  (a-c)  and 
time  series  of  lake  surface  elevation  (d). 


5.3.1  Materials  and  Methods 

Water  fluxes  to  and  from  Twelvemile  Lake  were  estimated  via  a  scoping  analysis  for  variable 
partitioning  of  precipitation  between  snow  and  rain,  different  flowpaths  of  snowmelt,  and 
historic  changes  in  growing  season  air  temperature  (May-September).  The  past  and  present 
spatial  distribution  of  permafrost,  which  acts  as  an  aquiclude  at  the  spatial  scale  necessary  for 
this  study,  is  not  completely  known;  therefore,  a  variety  of  permafrost  configurations  and 
associated  flow  pathways  are  considered  in  order  to  address  different  plausible  boundaries  of  the 
hydrologic  system.  Some  water  fluxes  to/from  the  lake  occur  directly  from  the  atmosphere 
(direct  flux),  while  other  fluxes  occur  indirectly  as  groundwater  or  overland  flow  (indirect  flux). 
Direct  fluxes  are  defined  to  be  snow,  rainfall,  and  lake  evaporation.  Indirect  fluxes  (Q,)  are 
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enumerated  as  follows  (Figure  5. 3. 1.1):  (i)  flows  from  the  contributing-area  slopes  of  the 

3  1 

watershed  {Qcomrib,  L  T'  ),  (ii)  subsurface  flow  across  the  watershed  boundary  over  a  “low  spot” 

3  1 

(“threshold”)  formed  by  a  topographically  depressed  permafrost  table  (Qthresh,  L  T'  ),  and  (ill) 
flows  to  or  from  the  lake  through  an  open  talik  (Qtaiik,  L  T'  ).  Indirect  fluxes  are  determined 
based  on  the  hydraulic  forces  and  hydraulic  properties  along  multiple  assumed  flowpaths.  The 
spatial  distribution  of  land  features  and  water  surface  elevations  are  determined  using  aerial 
photos  and  satellite  images  and  a  LIDAR  2.5  m-resolution  DEM  (Gesch,  2007).  Constraints  on 
the  current  spatial  distribution  of  permafrost  and  water  table  elevation  are  determined  from 
angering  and  frost  probing  during  2010  and  2011.  Mapped  resistivity  from  an  airborne 
electromagnetic  (AEM)  survey  conducted  during  2010  helped  to  visualize  the  overall  pattern  of 
ground  ice  in  the  watershed. 


(a)  View  looking  northeast 

Threshold, 
outlet  channel 


P.  SWE  E 

i  t 


Threshold, 
inlet  channel 


Ground 

surface 


(b)  View  looking  northwest 


P.  SWE  E  Ground  surface 


/Permafrost 


Figure  5.3. 1,1  Conceptual  model 
of  water  fluxes  considered  in  this 
study. 

View  looking  northeast, 
perpendicular  to  the  inlet  and 
outlet  channels  (a),  and  view 
looking  northwest  (h). 


Time  periods  of  1950-1980  and  1980-2010  were  considered,  the  former  period  assumed  to 
represent  conditions  prior  to  onset  of  lake  lowering  beginning  in  the  early  1980’s.  Years  1984 
and  2010  were  used  to  represent  lake  geometries  during  the  earlier  and  latter  periods, 
respectively.  Here,  an  indirect  flux  was  defined  to  be  significant  if  it  alone  provides  at  least  half 
of  the  indirect  flux  required  to  produce  the  observed  average  rate  of  lake  level  change  between 
1950  and  2010.  The  lake  water  budget  was  considered  to  have  significant  sensitivity  to  change  in 
a  particular  water  flux  if  that  change  would  produce  at  least  half  of  the  observed  change  in  lake 
lowering  beginning  in  the  early  1980’s. 

5.3.2  Results  and  Discussion 

Permafrost  was  not  found  via  frost  probing  and  soil  pit  analysis  at  19/21  sites  within  the  1984 
lake  boundary,  and  was  found  at  0.4-2. 5  m  depth  at  25/36  sites  outside  of  the  1984  lake  boundary 
(Figure  5.3.2. 1  a).  Most  of  the  ice-free  sites  outside  of  the  1984  lake  boundary  are  in  the  inlet  and 
outlet  channels,  where  ground  penetration  near  the  watershed  boundary  was  limited  due  to  the 
presence  of  shallow  gravel  and  water  (Jepsen  et  ah,  2012).  These  observations  indicate  that  the 
early  1980’s  lake  margin  approximately  separates  watershed  areas  with  and  without  permafrost. 
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excluding  the  channels.  This  pattern  is  consistent  with  that  from  the  AEM  resistivity 
interpretation  (Figure  5.3.2. 1  b),  where  resistivities  greater  than  about  530  ohm-m  (i.e.,  orange 
color  in  Figure  5.3.2. 1  b)  generally  reflect  the  presence  of  permafrost.  Based  on  these 
observations,  the  outer  and  inner  contributing  areas  used  for  the  Qcontrih  calculation  was  defined 
to  be  the  areas  outside  and  inside,  respectively,  the  1984  lake  boundary. 
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Figure  5,3,2, 1  Maps  showing  resistivity  of  subsurface  (upper  10  m)  from  the  AEM  survey 
and  depth  to  ice  from  ground-based  observations  (a),  and  AEM-based  resistivities  and 
lithology-ice  interpretations  along  cross  section  AA’  through  Twelvemile  Lake  (b).  Frozen 
material  in  (b)  indicates  permafrost.  Vertical  exaggeration  2,6x  in  panel  (b).  The  boundary 
separating  the  outer  (A  =  0)  and  inner  (A  =  1)  contributing  areas  is  shown  in  (a). 

The  increase  in  lake  evaporation  since  the  early  1980’s  due  to  increased  May-September  air 
temperatures  was  calculated  to  be  2  cm  yr'\  which  is  an  order  of  magnitude  lower  than  the 
observed  17.5  cm  yr'^  change  in  mass  balance,  indicating  that  increased  air  temperature  warming 
alone  is  not  likely  to  be  a  primary  mechanism.  Decreased  duration  of  the  ice-covered  season  may 
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be  as  much  as  1-2  weeks  since  the  late  1970’s  due  to  springtime  warming  (Brown  et  ah,  2010; 
Foster  et  ah,  2008).  This  would  increase  the  duration  of  the  open-water  season,  and  hence  lake 
evaporation  (Woo,  1980),  by  perhaps  10%,  which  is  still  not  adequate  to  be  a  primary 
mechanism. 

Analytical  results  using  the  equations  and  parameters  provided  in  detail  in  Jepsen  et  al.  (2013a) 
demonstrate  that  the  following  mechanisms  are  possible  explanations  for  the  lowering  of 
Twelvemile  Lake  beginning  in  the  early  1980’s:  (1)  changes  in  lateral,  shallow  (<  -40  m) 
groundwater  flow  across  the  watershed  boundary,  through  gravel,  (2)  changes  in  groundwater 
flow  resulting  from  the  development  of  an  open  talik,  provided  that  sub-permafrost  flow  is 
through  sand,  layers  of  which  are  sparse  relative  to  lacustrine  silt,  (3)  a  decrease  in  the  mass  of 
maximum-accumulation  snowpacks  in  conjunction  with  increased  infiltration  of  snowmelt  into, 
and  subsequent  evaporation  from,  fine-grained  sediment  around  the  lake. 

Increased  lateral,  shallow  groundwater  flow  from  a  lake’s  watershed  (mechanism  1)  would  occur 
if  a  lake  overtopped  the  permafrost-table  threshold  at  the  watershed  boundary,  resulting  in  a  type 
of  shallow-subsurface  flooding  (Michel  and  van  Everdingen,  1994).  This  process  would 
represent  a  subsurface  analogue  of  lake  drainage  via  overland  flow  through  surface  channels,  as 
observed  in  other  studies  (Woo  and  Mielko,  2007;  Brewer  et  ah,  1993;  Mackay,  1992),  and  may 
lead  to  sustained  lake  lowering  if  the  rate  of  thermal  erosion  at  the  permafrost  threshold  could 
equal  or  exceed  the  rate  of  lake  lowering  (Marsh  and  Neumann,  2001).  Thawing  of  the  transient 
layer  (permafrost  that  freezes  and  thaws  on  decadal  time  scales)  (Shur  et  ah,  2005)  in  gravelly 
material  could  facilitate  onset  of  this  process.  This  mechanism  was  found  to  be  potentially 
important  despite  the  low  topographic  gradients  of  Yukon  Flats  because  of  the  high  hydraulic 
conductivity  of  gravel. 

The  requirement  of  high  permeability  in  the  subpermafrost  aquifer  for  open-talik  development 
(mechanism  2)  to  be  an  important  mechanism  of  lake  lowering  results  from  the  low  hydraulic 
gradients  in  Yukon  Flats,  which  are  2-3  orders  of  magnitude  lower  than  gradients  observed  near 
Council,  Alaska,  where  open  taliks  are  an  important  mechanism  for  thermokarst  pond  drainage 
(Yoshikawa  and  Hinzman,  2003).  In  the  Yukon  Flats,  subpermafrost  sand  layers  are  minor  in 
abundance  relative  to  lacustrine  silt  (Williams,  1962;  Clark  et  ah,  2009),  suggesting  that  lake 
drainage  through  open  taliks  on  a  large  scale  in  Yukon  Flats  may  be  limited  by  the  permeability 
of  the  subpermafrost  aquifer. 

In  regard  to  increased  snowmelt  infiltration  (mechanism  3),  previous  studies  have  revealed  a 
wide  range  of  infdtration  capacities  of  frozen  soil  (for  reviews,  see  Ford  and  Bedford,  1987; 
Dingman,  1975).  In  environments  where  frozen  soil  has  been  found  to  prevent  snowmelt 
infiltration  and  enhance  overland  flow  (Woo  and  Steer,  1983;  Granger  et  ah,  1984;  Kane  and 
Stein,  1983;  Kane,  1980),  which  may  be  associated  with  relatively  cold,  moist  soils,  lake  mass 
would  have  a  high  sensitivity  to  snow  water  equivalent  (SWF).  In  contrasting  environments 
where  much  or  all  snowmelt  can  infiltrate  into  frozen  soil  (Suzuki  et  ah,  2006;  Carey  and  Woo, 
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1999),  lake  mass  would  have  low  sensitivity  to  SWE  beeause  of  how  little  snowmelt  reaches  the 
lake.  This  would  be  especially  true  in  areas  of  gentle  topography  and  fine-grained  soils,  such  as 
Yukon  Flats.  Progressive  lake  lowering  produces  a  gradient  in  the  exposed  land  surface  in 
vegetation  and  permafrost  distribution.  How  this  process  affects  soil  moisture  and  snowmelt 
infiltration,  and  possible  feedbacks  on  lake  recharge,  merits  further  study. 

All  of  the  plausible  mechanisms  for  lake  level  lowering  identified  in  this  study  require  the 
presence  of  certain  geologic  and/or  permafrost  conditions,  known  to  display  high  spatial 
variability  across  the  Yukon  Flats.  This  finding  is  consistent  with  the  high  degree  of  non¬ 
uniformity  observed  in  lake  level  changes  in  Yukon  Flats  (Rover  et  ah,  2012;  Riordan  et  ah, 
2006).  Examples  of  these  spatially  varying  characteristics  include  discontinuity  of  permafrost, 
depth  to  the  permafrost  table,  and  SWE  at  maximum  accumulation  (Sturm  et  ah,  2010,  2011). 


5.4  Lake  Talik  Evolution 

In  cold  regions,  hydrologic  systems  possess  seasonal  and  perennial  ice-free  zones  (taliks)  within 
areas  of  permafrost  that  control  and  are  enhanced  by  groundwater  flow  as  demonstrated  by 
model  simulations  described  in  sections  5.1  and  5.2.  To  further  develop  that  concept  within  the 
framework  of  our  study  area  on  the  lake  watershed  scale,  a  series  of  simulations  were  designed 
to  understand  talik  development  that  follows  lake  formation  in  watersheds  modeled  after 
conditions  in  the  Yukon  Flats.  The  goal  was  to  provide  insight  on  the  coupled  interaction 
between  groundwater  flow  and  ice  distribution.  The  modified  SUTRA  code  with  freeze-thaw 
physics  was  used  to  examine  timescales  and  controlling  factors  of  talik  formation.  This  work  was 
the  first  numerical  study  published  (Wellman  et  ah,  2013)  that  addressed  both  sub-lake  talik 
evolution  and  supra-permafrost  (active  layer)  near-lake  talik  evolution  under  the  influence  of 
seasonal  temperature  fluctuations  over  the  millennium  (1000  years)  timescale  for  a  range  of 
environmental  conditions.  Controlling  variables  explored  in  this  analysis  include  lake  depth  (and 
related  lake  size),  hydrologic  gradient  with  respect  to  the  lake,  and  shifts  in  climate-driven 
surface  temperature.  Information  derived  from  this  study  provides  insight  into  the  impact  of 
these  controls  on  expected  times  of  lake-talik  development,  on  shallow  and  deep  water  fluxes  to 
and  from  lakes,  and  on  lake-basin  permafrost  distribution,  thus  supplying  a  supporting 
foundation  for  important  scientific  linkages  such  as  biogeochemical  and  ecological  studies  of 
arctic  lake  changes  (Zimov  et  ah,  2006;  McGuire  et  ah,  2009). 

5.4.1  Materials  and  Methods 

The  basic  scenario  considered  in  this  modified-SUTRA  modeling  study  was  the  talik  evolution 
associated  with  a  newly-formed  lake  above  thick  (90-m)  permafrost.  The  model  enabled  thawing 
by  a  combination  of  thermal  conduction  and  advective  heat  transport.  Building  from  the  general 
climatic  and  structural  features  of  Twelvemile  Fake  (Figure  5.4.1),  a  suite  of  36  simulation 
scenarios  was  generated  to  address  variations  in  climate,  lake  size,  and  groundwater  flow 
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direction  on  talik  evolution.  The  suite  of  scenarios  allowed  examination  of  hydrologic  and 
cryologic  processes  under  a  range  of  conditions  that  typify  lakes  in  Yukon  Flats  and  in  other 
circumpolar  regions  of  relatively  thick  discontinuous  permafrost  and  low  topographic  relief. 


Figure  5,4,1  Cross-section  of  two-dimensional  cylindrical  model  with  boundary  conditions 
of  temperature,  fluid  pressure,  and  heat  flux  (a). 

Top  surface  consists  of  a  thermal  boundary  layer  between  air  temperature  (or  lake 
temperature)  and  aquifer  surface  temperature,  with  a  specified  air  (or  lake)  temperature 
on  the  top  surface  and  a  specified  pressure  on  the  lower  surface  representing  the  land 
surface  pressure  (p=0)  or  lake  bottom  pressure  (p=hydrostatic  pressure  at  local  lake 
depth).  Annual  synthetic  lake  bottom  temperature  as  a  function  of  various  lake  depths 
patterned  after  Burn  (2002,  2005)  (b). 
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Generalized  geologic  layer  representation  for  the  2-D  cylindrical  model  followed  the  basic 
lithology  of  the  Fort  Yukon  area  as  reported  in  Clark  et  al.  (2009)  and  Williams  (1962)  and  as 
interpreted  from  the  AEM  survey  of  the  region  (Minsley  et  ah,  2012).  A  thin  (2  m)  surficial 
sedimentary  deposit  characterized  as  silt  and  sand  exists  above  a  thick  (28  m)  gravel  and  sand 
layer,  which  is  underlain  by  an  extensive  undifferentiated  deposit  believed  to  be  predominately 
silt  in  composition. 

Constant  and  time -varying  conditions  that  specify  temperature,  fluid  pressure,  and  subsurface 
geothermal  energy  flux  were  applied  to  the  model  boundaries  (Figure  5.4.1).  The  ground  surface 
at  0  m  depth  consisted  of  an  idealized  lake  bottom  with  prescribed  fluid  pressures  and  oscillating 
seasonal  temperatures  reflective  of  an  overlying  lake  body 

Prescribed  head  (pressure)  conditions  in  the  sub-permafrost  aquifer  were  used  to  control  the 
direction  of  water  movement  between  the  sub-permafrost  groundwater  and  simulated  lake 
bottom  surface.  The  examined  scenarios  include:  (1)  hydrostatic  conditions  without  an  exchange 
of  water  between  the  lake  and  underlying  sub-permafrost  aquifer,  (2)  gaining  lake  conditions 
with  the  potential  for  sub-permafrost  groundwater  to  discharge  to  the  lake,  and  (3)  losing  lake 
conditions  with  the  potential  for  water  to  discharge  from  the  lake. 

Three  climates  were  examined,  reflective  of  the  present,  colder  than  present,  and  warmer  than 
present.  The  yearly  temperature  cycles  for  the  “colder”  and  “warmer”  climate  scenarios  were 
modified  from  that  of  the  current  climate  by  reducing  and  raising,  respectively,  temperatures  by 
2°C  for  the  terrestrial  areas  and  by  1°C  along  the  lake  bottom.  The  adjustments  of  temperature 
assumed  an  incremental  change  to  the  mean  conditions  while  the  amplitude  of  seasonal 
fluctuations  was  preserved.  The  rationale  for  exploring  different  climates  was  based  on 
documented  changes  in  historic  air  temperature  (Serreze  et  ah,  2000)  and  projected  GCM  (global 
climate  model)  changes  in  future  air  temperatures  in  Alaska  (Chapman  and  Walsh,  2007;  Walsh 
et  ah,  2008).  The  selected  2°C  temperature  adjustment  was  based  on  the  50-year  projected  GCM 
annual  average  temperature  increase  for  the  Yukon  Flats  assuming  mid-range  emissions 
(Scenarios  Network  for  Alaska  and  Arctic  Planning  (SNAP),  2012).  For  consistency,  an 
equivalent  reduction  in  temperature  for  a  climate  colder  than  present  was  also  explored 
indicative  of  conditions  half  a  century  or  longer  in  the  past. 

5.4.2  Results  and  Discussion 

A  detailed  comparison  of  transient  results  for  the  12-m  deep  losing,  gaining,  and  static  lake 
conditions,  all  under  the  current  climate,  provides  an  opportunity  for  describing  the  inter¬ 
relations  and  temporal  patterns  of  ice  structure  and  groundwater  flow.  Figure  5.4.2. 1  a,b,  and  c 
illustrates  the  temporal  progression  in  flow  and  ice  volume  through  time  since  lake  formation  for 
the  losing,  gaining,  and  static  lake  conditions,  respectively.  The  amplitude  of  yearly  fluctuations 
in  ice  volume  and  groundwater  flow  are  expressed  in  Figure  5.4.2. 1  as  the  vertical  thickness  of 
the  color  bands,  because  the  individual  yearly  fluctuations  appear  as  a  continuous  band  on  the 
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long  time  scale  plotted.  Ice  volume  present  under  the  lake  decreases  with  time.  Sub-lake  talik 
breakthrough  is  coincident  with  a  break  in  the  slope  indicating  the  change  in  ice  volume  with 
time  for  regions  with  >0.1  ice  saturation  and  is  denoted  by  the  vertical  dashed  line  in  Figure 
5.4.2. 1.  After  breakthrough,  the  sub-lake  talik  expands  laterally  outward  from  the  lake,  further 
reducing  the  ice  volume  under  the  lake  through  time.  Sub-lake  talik  breakthrough  occurs  in  215 
years  for  the  losing  lake,  240  years  for  the  gaining  lake,  and  598  years  for  the  static  lake, 
yielding  an  effective  rate  of  sub-lake  vertical  thaw  of  approximately  0.42  m  yr‘\  0.38  m  yr'\  and 
0.15  m  yr'^  respectively.  The  rate  of  total  ice  volume  decrease  (approximated  by  the  decline  in 
ice  volume  with  an  ice  saturation  >  0.10),  and  the  sub-lake  groundwater  flow  both  increase  once 
an  open  talik  is  developed  for  the  losing  and  gaining  systems,  providing  a  feedback  mechanism 
that  enhances  the  thaw  rate  via  heat  advection.  Prior  to  this  point,  reduction  in  ice  volume  is 
largely  controlled  by  heat  conduction,  because  groundwater  flow  is  restricted  by  pore  ice  that 
impedes  groundwater  flow.  Yearly  fluctuations  in  ice  volume  diminish  with  time  for  the  losing 
and  gaining  lakes  thereafter  as  the  ice  structure  along  the  lake  margin  tends  toward  equilibrium. 

Four  cases  of  talik  evolution  are  illustrated  for  a  12-m  lake:  hydrostatic  and  gaining  lakes  under 
the  current  climate  (Figure  5. 4. 2.2  a,  b),  and  gaining  and  losing  lakes  under  a  warmer  climate 
(Figure  5. 4.2. 2  c,  d).  Results  are  examined  at  100,  250,  and  1000  year  elapsed  times  representing 
early  to  late  phases  in  talik  development.  At  100  years,  thaw  in  all  of  the  scenarios  examined  is 
conduction-controlled,  because  ice  saturation  is  sufficient  to  impede  groundwater  flow 
irrespective  of  whether  a  hydraulic  gradient  exists  between  the  lake  and  aquifer  or  if  a  warmer 
climate  is  imposed.  At  250  years  elapsed  time,  the  hydrostatic  lake  has  a  thaw  bulb  that  extends 
to  a  depth  of  about  40  m  at  the  lake  center.  For  a  gaining  lake  with  current  and  warmer  climates. 
Figure  5. 4. 2.2  b  and  c,  respectively,  permafrost  thaw  is  enhanced  relative  to  the  hydrostatic  lake 
situation  due  to  advection  of  heat  in  the  upward-flowing  groundwater.  The  situation  with  a  losing 
lake  under  the  warmer  climate,  exhibits  the  fastest  thawing  among  the  four  12  m  lake  cases 
considered  (Figure  5. 4. 2.2  d).  The  talik  structure  is  advanced  within  250  years  and  possesses  a 
width  nearly  four  times  greater  than  that  of  the  gaining  lake  under  the  current  climate.  After  1000 
years,  the  system  exhibits  extensive  permafrost  loss  with  the  thaw  front  extending  200  m  beyond 
the  lake  margin.  The  reduction  in  permafrost  within  the  lake  watershed  is  due  to  warmer  surface 
temperatures  that  couple  with  flowing  groundwater  to  alter  both  the  temperature  distribution  and 
ice  saturation. 
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—  -  Sub-permafrost  How 

-  Ice  under  lake  (Ice  sat.  >  0.9) 

-  Ice  under  lake  (Ice  sat.  >0.1) 


Figure  5,4,2, 1  Temporal  evolution  of  groundwater  flow  and  ice  volume  for  the  case  of  a 
12m  lake  that  is  losing  (a),  gaining  (h),  and  hydrostatic  (c)  for  the  current  climate. 

The  left  axis  measures  groundwater  flow  (>  0  groundwater  enters  aquifer,  <  0  groundwater 
leaves  aquifer)  and  the  right  axis  measures  the  volume  of  ice  for  ice  saturations  >0,1  and 
>0,9,  The  dashed  line  in  a-c  indicates  the  talik  breakthrough  time  for  each  case.  Broad 
color  hands  result  from  yearly  oscillations  of  the  plotted  factor. 
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Figure  5,4,2,2  Cross-sections  of  2D  cylindrical  model  (mirrored)  showing  states  of  lake- 
talik  formation,  ice  saturation,  and  groundwater  movement  where  a  12m  lake  forms  over  a 
continuous  permafrost  layer  of  about  90m  in  thickness. 

Figure  columns  distinguish  early,  middle,  and  late-time  conditions  at  100,  250,  and  1000 
years,  respectively,  from  the  time  of  lake  formation.  Figure  rows  differentiate  four 
scenarios:  (a)  current  climate  with  hydrostatic  (no  flow)  lake,  (h)  current  climate  with 
gaining  lake,  (c)  warmer  climate  with  gaining  lake,  and  (d)  warmer  climate  with  losing 
lake. 

Average  times  of  sub-lake  talik  breakthrough  for  the  different  climates,  directions  of  lake  flow, 
and  lake  depths  reveal  relational  patterns  (Figure  5.4. 2. 3  a).  For  processes  exceeding  the  time  of 
analysis,  the  maximum  analysis  time  is  used  to  determine  group  averages,  thereby  yielding 
minimum  estimates.  For  the  warmer  climate  simulations,  talik  breakthrough  times  are  shorter 
than  for  cooler  climates.  Time  required  for  open  talik  formation  decreases  by  -30  %  from  the 
colder  climate  scenarios  to  the  warmer  climate  scenarios  on  average.  Talik  breakthrough  time 
can  vary  by  more  than  a  factor  of  five  between  hydrostatic  (requiring  >1000  years  for  some 
cases)  and  gaining  or  losing  conditions  (forming  as  quickly  as  -200  years)  due  to  the  additional 
component  of  heat  advection  for  the  latter  cases.  The  presence  or  absence  of  groundwater  flow 
to/from  the  sub-permafrost  aquifer  also  exerts  a  large  influence  on  equilibrium  times  of 
lake/groundwater  exchange  (Figure  5. 4. 2. 3  b).  In  the  case  of  a  hydrostatic  lake,  inflow  is  derived 
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only  from  the  supra-permafrost  aquifer,  and  equilibrium,  sensitive  to  proeesses  near  the  lake 
margin,  is  quiekly  aehieved  in  just  1  to  9  years  (<  deeade).  Times  of  flow  equilibrium  for  losing 
and  gaining  lakes  are  mueh  longer  and  average  661  and  694  years,  respeetively. 

Caleulated  times  to  aehieve  equilibrium  in  permafrost  volume  average  606  and  711  years  for 
permafrost  lee  saturations  of  0.90  and  0.10,  respeetively,  with  a  range  between  206  to  >  900 
years.  In  terms  of  group  averages,  elimate  has  minimal  impaet  on  eryologie  equilibrium  times  for 
regions  of  high  lee  saturation  permafrost  (lee  saturation  >  0.9),  and  a  moderate  impaet  for  the 
bulk  of  the  permafrost  volume  (lee  saturation  >  0.1)  (Figure  5.4. 2. 3  e,  d). 

Comparison  between  groundwater  flux  within  the  supra-permafrost  aquifer  and  net  flux  of  water 
to  the  lake  at  equilibrium  exhibits  stark  differenees  depending  on  system  eonditions  (Figure 
5. 4. 2. 3  e,  f).  Groundwater  exehange  with  the  lake  through  the  open  talik  is  greater  than  flow 
within  the  supra-permafrost  aquifer  by  more  than  -'1.5  orders  of  magnitude,  on  average.  For  the 
gaining  or  losing  lakes  evaluated,  flow  through  the  supra-permafrost  aquifer  is  relatively  small 
(typieally  <  1  %)  eompared  to  flow  exehanged  between  the  lake  and  sub-permafrost  aquifer. 
However,  for  a  hydrostatie  lake,  supra-permafrost  flow  is  the  dominant  groundwater  eontribution 
to  the  lake. 

Climate  warming  eauses  flow  within  the  supra-permafrost  aquifer  to  inerease  by  approximately  2 
orders  of  magnitude  from  10  to  10  m  yr'  aeross  the  speetrum  of  elimate  regimes  eonsidered 
(Figure  5. 4. 2. 3  f).  This  trend  results  from  thiekening  of  the  supra-permafrost  aquifer  and  thus 
greater  groundwater  exehange  aeross  a  more  developed  lake  margin  that  is  permeable  for  a 
prolonged  part  of  the  year  for  progressively  warmer  seenarios. 

The  results  demonstrate  the  influenee  of  time,  elimate,  lake  size,  and  hydraulie  eonditions  on 
talik  development  for  eold-region  lakes.  Regardless  of  the  eondition  evaluated,  early-time  thaw 
after  the  formation  of  a  lake  (periods  of  100  years  or  less)  is  eontrolled  by  heat  eonduetion  below 
the  lake  bed.  At  later  times,  the  transfer  of  heat  by  groundwater  flow,  where  signifieant 
groundwater  flow  oeeurs,  is  an  important  meehanism  that  aeeelerates  permafrost  thaw  and 
expedites  breakthrough  of  sub-lake  taliks.  Warmer  elimates  generally  inerease  thaw  rates  and  for 
the  systems  examined  exhibit:  (a)  thiekening  of  and  inereased  flow  through  the  supra-permafrost 
aquifer,  (b)  reduetion  of  sub-lake  open  talik  formation  times  and  more  rapid  talik  expansion, 
thereby  enhaneing  earlier  hydraulie  eommunieation  between  the  lake  and  sub-permafrost  aquifer, 
and  (e)  greater  overall  loss  of  permafrost. 

Hydrologie  eonneetions  between  the  lake  and  sub-permafrost  aquifer  that  depart  from 
hydrostatie  eonditions  are  expeeted  to  allow  intense  transfers  of  water  and  energy  for  short  time 
periods  in  oases  suoh  as  sub-lake  talik  formation.  Given  that  lakes  in  Yukon  Flats  are  observed  to 
ohange  in  size,  appear,  disappear,  and  even  reappear  in  pre-existing  looations,  brief  periods  of 
thawing,  refreezing  and  oooasional  intense  vertioal  groundwater  exehange  may  be  a  reourring 
prooess  throughout  the  region. 
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Figure  5,4,2.3  Group  averages  of  simulation  results  for  three  categories:  climate  regime, 
lake  flow,  and  maximum  lake  depth. 

Tabulated  results  for  the  36  simulations  are  presented  in  Table  3  of  Wellman  et  al.,  2013. 
Results  are  shown  for:  talik  breakthrough  time  (a),  time  to  reach  flow  equilibrium  at  the 
lake  bottom(b),  time  to  reach  equilibrium  in  residual  permafrost  volume  for  ice  saturation 
>0.9  and  >0.1,  respectively(c,  d),  equilibrium  groundwater  flow  across  the  lake  bottom  (e), 
and  mean  annual  groundwater  flow  to  the  lake  through  supra-permafrost  aquifer  (f). 
Grey-shaded  bars  show  averages  for  each  analysis  and  thin  black  range  bars  indicate 
minimum  and  maximum  values  among  all  analyses. 
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In  summary,  although  it  is  generally  assumed  that  climate  is  the  dominant  environmental  driver 
in  cold  regions,  lake  size  and  groundwater  flow  are  shown  to  impose  equally  important  controls 
on  evolution  of  permafrost  thaw,  timing  of  talik  formation,  and  lake/groundwater  exchange.  The 
distribution  of  permafrost  within  a  cold-region  lake  watershed  will  influence,  and  be  influenced 
by,  lake-groundwater  exchange  through  supra-permafrost  and  sub-lake  taliks.  Studies  that 
evaluate  changes  in  lake  size  and  distribution  and  that  predict  climate-change-driven  hydrologic 
evolution  in  cold  regions  should  consider  groundwater-ice  dynamics. 


5,5  Sensitivity  Analysis  of  Geophysical  Data  to  Permafrost  Conditions 

The  general  hydrogeologic  framework  for  the  coupled  thermal-hydrologic  simulations,  described 
above,  was  based  largely  on  interpretations  from  electromagnetic  geophysical  information  in  the 
Yukon  Flats.  However,  some  of  the  geophysics-derived  interpretations,  including  lake-talik 
conditions,  have  very  limited  ground  truth  information  for  calibration  and  validation  due  to  the 
paucity  of  standard  hydrogeologic  data  collected  in  this  remote  region.  To  address  potential 
sources  of  uncertainty  associated  with  geophysical  interpretations,  we  investigated  the  ability  of 
geophysical  measurements  to  recover  information  about  the  underlying  spatial  distribution  of 
permafrost  and  hydrologic  properties.  This  was  accomplished  in  three  steps:  (1)  development  of 
a  physical  property  relation  that  connects  permafrost  and  hydrologic  properties  to  geophysical 
properties;  (2)  simulation  of  synthetic  geophysical  data  that  would  be  expected  for  various 
permafrost  hydrologic  conditions;  and  (3)  inversion  of  the  synthetic  geophysical  data  using 
realistic  levels  of  noise  to  investigate  the  ability  to  resolve  specific  physical  features  of  interest. 
The  focus  was  on  electromagnetic  geophysical  methods  as  these  types  of  data  had  previously 
been  acquired  near  Twelvemile  Lake  in  the  Yukon  Flats,  Alaska  (Ball  et  ak,  2011;  Minsley  et  al., 
2012).  Twelvemile  Lake  is  the  lake  that  was  the  basis  for  the  lake  simulations  described  in  the 
previous  section. 

5.5.1  Materials  and  Methods 

Model  output  from  the  lake  talik  evolution  simulations  (section  5.4)  utilizing  the  modified 
SUTRA  code  was  used  as  the  basis  for  generating  synthetic  geophysical  data.  Relevant  SUTRA 
output  including  temperature,  pressure,  and  ice  saturation  at  each  simulation  time  step  for  each 
model  scenario  was  converted  to  electrical  resistivity,  the  geophysical  property  needed  to 
simulate  electromagnetic  data,  as  follows. 
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Bulk  electrical  conductivity  is  described  by  Revil  (2012)  as 


(j 


F 


af  +  m(s-J'F-\)a^ 


where  a  is  the  bulk  electrical  conductivity  [S/m];  is  the  fractional  water  saturation  [-]  in  the 
pore  space,  where  Sw=  1-  S,  and  S,  is  the  fractional  ice  saturation  [-]  in  the  pore  space;  Of  is  the 
conductivity  of  the  saturating  pore  fluid  [S/m];  m  is  the  Archie  cementation  exponent  [-];  n  is  the 
Archie  saturation  exponent  [-];  F  is  the  formation  factor  [-],  where  F  =  and  (j)  is  the  matrix 
porosity  [-];  and  is  the  conductivity  [S/m]  associated  with  grain  surfaces.  The  Archie 
exponents  m  and  n  are  known  to  vary  as  a  function  of  pore  geometry;  here,  we  use  m  =  n  =  1.5. 
Simulation  results  are  presented  as  electrical  resistivity  [ohm-m],  which  is  the  inverse  of  the 
conductivity,  i.e.  p  =  1/a. 


Information  about  the  different  lithologic  units  incorporated  in  the  model  analysis  (Figure  5.5.1) 
was  used  to  define  static  model  properties  such  as  porosity,  grain  mass  density,  cation  exchange 
capacity,  and  Archie’s  exponents  (Table  5.5.1).  Dynamic  outputs  from  the  SUTRA  simulations, 
including  temperature  and  ice  saturation,  were  combined  with  the  static  variables  to  predict  the 
evolving  electrical  resistivity  structure. 


Figure  5,5,1  Axis-symmetric  model  geometry  indicating  different  lithologic  units  and 
simulated  lake  depths/extents. 
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Geologic  unit  properties 

Lithology; 

Unit  1 

Sediment  (silty  sand) 

Unit  2 

Sediment  (gravelly  sand) 

Unit  3 

Lacustrine  silt 

Unit  depth  range  [m] 

0-2,  2-30,  30-250 

Porosity  [-] 

0.25,  0.25,  0.20 

Geophysical  parameters 

Archie  cementation  exponent  (m)  [-] 

1.5 

Archie  saturation  exponent  (n)  [-] 

1.5 

Water  salinity  (Q  [ppm] 

250  (Si  =  0) 

Na  ionic  mobility  (P)  [m  W s] 

5.8  X  10'^  (25°C) 

Cr  ionic  mobility  (P)  [m^Ws] 

7.9  X  10'^  (25°C) 

Na  surface  ionic  mobility  (Ps)  [m  W s] 

0.51  X  10'^  (25°C) 

Grain  mass  density  (pg)  [kg/m  ] 

2650 

Cation  exchange  capacity  (x)  [C/kg] 

200,  10,  500 

Table  5,5.1  Description  of  geologic  units  and  physical  properties  used  in  numerical 
simulations.  Entries  separated  by  commas  represent  parameters  with  different  values  for 
each  of  the  lithologic  units. 

Synthetic  airborne  electromagnetic  (AEM)  data  were  simulated  for  each  snapshot  of  predicted 
bulk  resistivity  values  using  nominal  system  parameters  based  on  the  Fugro  RESOLVE 
frequency-domain  AEM  system  that  was  used  in  the  Yukon  Flats  survey  (Minsley  et  ah,  2012). 
Data  are  simulated  at  the  nominal  survey  elevation  of  30  m  above  ground  surface  using  the  one¬ 
dimensional  modeling  equations  described  in  Minsley  (2011).  A  200-m  vertical  profile  of 
resistivity  as  a  function  of  depth  was  extracted  at  each  survey  location  and  was  used  to  simulate 
forward  geophysical  responses.  Forward  simulations  were  repeated  for  each  of  the  50  simulation 
times  between  0  and  1,000  years  output  from  SEITRA,  resulting  in  9,050  data  locations  per 
modeling  scenario. 

Synthetic  ground-based  electromagnetic  data  were  simulated  using  nominal  system  parameters 
based  on  the  GEM-2  instrument  (Huang  and  Won,  2003).  A  system  elevation  of  1  m  above 
ground  was  assumed,  which  is  typical  for  this  hand-carried  instrument. 
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To  address  the  inherent  uncertainty  associated  with  geophysical  inversion,  a  Bayesian  Markov 
chain  Monte  Carlo  (McMC)  algorithm  developed  for  frequency-domain  EM  data  (Minsley, 
2011)  was  utilized  to  explore  the  ability  of  simulated  AEM  data  to  recover  the  true  distribution 
of  subsurface  resistivity  values  throughout  the  1,000-year  lake  talik  simulations.  At  every  data 
location  along  the  survey  profile,  an  ensemble  of  100,000  resistivity  models  was  sampled 
according  to  the  Metropolis-Hastings  algorithm  (Hastings,  1970;  Metropolis  et  ah,  1953). 
According  to  Bayes’  theorem,  each  model  was  assigned  a  posterior  probability  that  was  a 
measure  of  (1)  its  prior  probability  which,  in  this  case,  was  used  to  penalize  models  with 
unrealistically  large  contrasts  in  resistivity  over  thin  layers;  and  (2)  its  data  likelihood,  which  was 
a  measure  of  how  well  the  predicted  data  for  a  given  resistivity  model  match  the  observed  data 
within  data  errors.  Numerous  measures  and  statistics  were  generated  from  the  ensemble  of 
plausible  resistivity  models,  such  as:  the  single  most-probable  model,  the  probability  distribution 
of  resistivity  values  at  any  depth,  the  probability  distribution  of  where  layer  interfaces  occur  as  a 
function  of  depth,  and  the  probability  distribution  of  the  number  of  layers  (model  complexity) 
needed  to  fit  the  measured  data.  Einally,  probability  distributions  of  resistivity  were  combined 
with  assumptions  about  the  distribution  of  resistivity  values  for  any  lithology  and/or  ice  content 
in  order  to  make  a  probabilistic  assessment  of  lithology  or  ice  content. 

5.5.2  Results  and  Discussion 

Eor  each  of  the  1,000-year  lake  talik  evolution  simulations,  the  static  variables  summarized  in 
Table  5.5.1  were  combined  with  the  spatially  and  temporally  variable  state  variables  T  and  Si 
output  by  SUTRA  to  predict  the  distribution  of  bulk  resistivity  at  each  time  step.  An  example  of 
SUTRA  output  variables  for  the  6  m-deep  gaining  lake  scenario  at  240  years  (the  approximate 
sub-lake  talik  breakthrough  time  for  that  scenario)  is  shown  in  Eigure  5.5.2. 1  A-B,  and  the 
predicted  resistivity  for  this  simulation  step  is  shown  in  Eigure  5.5.2. 1  C.  The  influence  of 
different  lithologic  units  is  clearly  manifested  in  the  predicted  resistivity  values,  whereas 
lithology  is  not  overly  evident  in  the  SUTRA  state  variables.  Eor  a  single  unit,  there  is  a  clear 
difference  in  resistivity  for  frozen  versus  unfrozen  conditions.  Across  different  units,  there  is  a 
contrast  in  resistivity  when  both  units  are  frozen  or  unfrozen.  Resistivity  can  therefore  be  a 
valuable  indicator  of  both  geologic  and  ice  content  variability.  However,  there  is  also  ambiguity 
in  resistivity  values  as  both  unfrozen  Unit  2  and  frozen  Unit  3  appear  to  have  intermediate 
resistivity  values  of  approximately  100-300  ohm-m  (Eigure  5.5.2. 1  C)  and  cannot  be 
characterized  by  their  resistivity  values  alone.  This  ambiguity  in  resistivity  can  only  be 
overcome  by  additional  information  such  as  borehole  data  or  prior  knowledge  of  geologic 
structure. 
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Figure  5,5,2, 1  SUTRA  model  outputs  and  geophysical  transformations  from  the  6-m 
gaining  lake  simulation  at  240  years.  Ice  saturation  (A)  and  temperature  (B)  are  converted 
to  predictions  of  hulk  resistivity  (C),  Variahility  in  resistivity  as  a  function  of  temperature 
is  indicated  in  (D)  for  lithologic  units  1-3, 

Lithology  and  ice  saturation  are  the  primary  factors  that  control  simulated  resistivity  values, 
though  ice  saturation  is  a  function  of  temperature.  The  empirical  relation  between  temperature 
and  bulk  resistivity  is  shown  in  Figure  5.5.2. 1  D  by  cross-plotting  values  from  Figure  5.5.2. 1  B- 
C.  Within  each  lithology,  resistivity  is  relatively  constant  above  zero  degrees,  with  a  rapid 
increase  in  resistivity  for  temperatures  below  zero  degrees.  This  result  is  very  similar  to  the 
temperature-resistivity  relationships  illustrated  by  Hoekstra  et  al.  (1975;  Fig.  1),  lending 
confidence  to  the  physical  property  definitions  used  in  this  study.  Above  zero  degrees,  the  slight 
decrease  in  resistivity  is  due  to  the  temperature-dependence  of  fluid  resistivity.  The  rapid 
increase  in  resistivity  below  zero  degrees  is  primarily  caused  by  reductions  in  effective  porosity 
due  to  increasing  ice  saturation,  though  changes  in  surface  conductivity  and  salinity  at  increasing 
ice  saturation  are  also  contributing  factors.  Below  -1  °C,  the  change  in  resistivity  values  as  a 
function  of  temperature  rapidly  decreases.  This  is  an  artifact  caused  by  the  imposed  temperature- 
ice  saturation  relationship  defined  in  SUTRA  that,  for  these  examples,  enforces  99%  ice 
saturation  at  -1  °C.  It  is  more  likely  that  ice  saturation  continues  to  increase  asymptotically  over 
a  larger  range  of  temperatures  below  zero  degrees,  with  corresponding  increases  in  electrical 
resistivity.  However,  because  AEM  methods  are  limited  in  their  ability  to  discern  differences 
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among  very  high  resistivity  values,  this  artifaet  does  not  significantly  impact  the  results 
presented  here. 

Geophysical  data  (not  shown)  were  simulated  for  each  of  the  electrical  resistivity  models.  The 
simulated  data  were  then  used  to  recover  estimates  of  the  original  resistivity  values.  For  example, 
resistivity  parameter  estimation  results  for  the  6-m  deep  hydrostatic  lake  scenario  are  shown  in 
Figure  5. 5. 2.2.  At  each  location  along  the  profile,  the  average  resistivity  model  as  function  of 
depth  was  calculated  from  the  McMC  ensemble  of  100,000  plausible  models.  The  overall  pattern 
of  different  lithologic  units  and  frozen/unfrozen  regions  was  accurately  depicted  with  two 
exceptions:  (1)  the  specific  distribution  of  partial  ice  saturation  beneath  the  lake  before  thaw  has 
equilibrated;  and  (2)  the  shallow  sand  layer  (Unit  1)  that  is  generally  too  thin  to  be  resolved 
using  AEM  data.  A  point-by-point  comparison  of  true  versus  predicted  resistivity  values  for  the 
hydrostatic  6  m-deep  lake  scenario  at  the  simulation  time  1,000  years  is  shown  in  Figure  5. 5. 2. 3. 
The  cross-plot  of  true  versus  estimated  resistivity  values  generally  fall  along  the  1:1  line, 
providing  a  more  quantitative  indication  of  the  ability  to  estimate  the  subsurface  resistivity 
structure.  Estimates  of  the  true  resistivity  values  for  each  lithology  and  freeze/thaw  state  (Eigure 
5. 5. 2. 3  B)  tend  to  be  indistinct;  appearing  as  a  vertical  range  of  possible  values  in  Eigure  5. 5. 2. 3 
A  due  to  the  inherent  resolution  limitations  of  inverse  methods  and  parameter  tradeoffs  (Day- 
Eewis,  et  ah,  2005;  Oldenborger  and  Routh,  2009).  Although  the  greatest  point  density  for  both 
frozen  and  unfrozen  silts  (Unit  3)  falls  along  the  1 : 1  line,  resistivity  values  for  these  components 
of  the  model  are  also  often  overestimated;  this  is  likely  due  to  uncertainties  in  the  location  of  the 
interface  between  the  silt  and  gravel  units.  This  is  in  contrast  with  the  systematic 
underestimation  of  frozen  gravel  resistivity  values  due  to  the  inability  to  discriminate  very  high 
resistivity  values  using  EM  methods  (Ward  and  Hohmann,  1988).  Erozen  sands  (true  log 
resistivity  -'2.8  in  Eigure  5. 5. 2. 3  B)  are  also  systematically  overestimated  in  Eigure  5. 5.2. 3  A;  in 
this  case,  due  to  the  inability  to  resolve  this  relatively  thin  resistive  layer. 
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Figure  5,5,2,2  Mean  resistivity  model  extracted  from  McMC  ensembles.  Results  are  shown 
for  the  6-m-deep  hydrostatic  lake  scenario  outputs  at  (A)  100  years,  (B)  240  years,  and  (C) 
1,000  years. 


While  useful,  single  ‘best’  estimates  of  resistivity  values  at  any  location  (Figure  5. 5. 2. 3)  are  not 
fully  representative  of  the  information  contained  in  the  AEM  data  and  associated  model 
uncertainty.  From  the  McMC  analysis  of  100,000  models  at  each  data  location,  estimates  of  the 
posterior  probability  density  function  (pdf)  of  resistivity  were  generated  for  each  point  in  the 
model.  Probability  distributions  were  extracted  from  a  depth  of  15  m,  within  the  gravel  layer 
(Unit  2),  at  one  location  where  unfrozen  conditions  exist  (r  =  0  m),  and  a  second  location  outside 
the  lake  extent  (r  =  750  m)  where  the  ground  remains  frozen  (Figure  5. 5. 2.4  A).  Results  from  a 
depth  of  50  m,  within  the  silt  layer  (Unit  3),  are  shown  in  Figure  5. 5.2. 4  B.  With  the  exception 
of  the  frozen  gravels,  whose  resistivity  tends  to  be  underestimated,  the  peak  of  each  pdf  is  a  good 
estimate  of  the  true  resistivity  value  at  that  location. 
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Figure  5.5.2,3  Performance  of  geophysical  parameter  estimation  in  recovering  true 
parameter  values.  True  versus  McMC-estimated  resistivity  values  for  the  hydrostatic  6-m- 
deep  lake  scenario  at  simulation  time  1,000  years  (A),  compared  with  the  frequency 
distribution  of  true  resistivity  values  (B).  Estimated  resistivity  values  generally  fall  along 
the  dashed  1:1  line  in  (A),  with  exceptions  being  under-prediction  of  the  resistive  frozen 
gravels,  over-prediction  of  the  thin  surficial  frozen  sand,  and  some  over-prediction  of  the 
frozen  silt  where  it  is  in  contact  with  frozen  gravel. 
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Figure  5,5.2,4  McMC-estimated  resistivity  posterior  distributions  within  frozen  and 
unfrozen  unit  #2  gravels  (A)  and  frozen  and  unfrozen  unit  #3  silts  (B)  for  the  hydrostatic  6- 
m-deep  lake  scenario  at  1,000  years.  Unfrozen  resistivity  distributions  are  extracted 
beneath  the  center  of  the  lake  (r  =  0)  at  depths  of  15  m  and  50  m  for  the  gravels  and  silts, 
respectively.  Frozen  distributions  are  extracted  at  the  same  depths,  but  at  r  =  750  m.  The 
upper  x-axes  labels  indicate  approximate  ice  saturation  based  on  the  lithology-dependent 
ice  saturation  versus  resistivity  curves. 

Further  illustration  of  the  spatial  and  temporal  ehanges  in  resistivity  pdfs  are  shown  in  Figure 
5. 5. 2. 5.  The  resistivity  pdf  is  displayed  as  a  funetion  of  distanee  from  the  lake  eenter  at  the  same 
depths  (15  m  and  50  m)  shown  in  Figure  5. 5. 2. 5,  eorresponding  to  gravel  (Figure  5. 5.2. 5  A,C, 
and  E)  and  silt  (Figure  5. 5. 2. 5  B,  D,  and  F)  loeations.  High  probabilities,  i.e.  the  peaks  in  Figure 
5. 5. 2.4,  eorrespond  to  dark-shaded  areas  in  Figure  5. 5. 2. 5.  Images  are  shown  for  three  different 
time  steps  in  the  SUTRA  simulation  for  the  hydrostatie  6  m-deep  lake  seenario:  100  years 
(Figure  5. 5. 2. 5  A-B),  240  years  (Figure  5. 5. 2. 5  C-D),  and  1,000  years  (Figure  5. 5. 2. 5  E-F). 

Approximate  iee-saturation  values,  translated  from  the  iee  versus  resistivity  relationships  for 
eaeh  lithology,  are  displayed  on  the  right  axis  of  eaeh  panel  in  Figure  5. 5.2. 5,  and  true  resistivity 
values  are  plotted  as  a  dashed  line.  Observations  from  Figure  5. 5. 2. 5  inelude: 

(1)  Outside  the  lake  boundary,  pdfs  are  signifieantly  more  sharply  peaked  (darker  shading)  for 
the  gravel  unit  than  the  silt  unit,  suggesting  better  resolution  of  shallower  resistivity  values 
within  the  gravel  layer.  It  should  be  noted  however,  that  this  improved  resolution  does  not  imply 
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improved  model  accuracy;  in  fact,  the  highest  probability  region  slightly  underestimates  the  true 
resistivity  value.  (2)  Probability  distributions  for  the  silt  layer  track  the  true  values,  but  with 
greater  uncertainty.  (3)  Inside  the  lake  boundary,  gravel  resistivity  values  are  not  as  well 
resolved  compared  with  locations  outside  the  lake  boundary  due  to  the  loss  of  signal  associated 
with  the  relatively  conductive  lake  water.  (4)  Increasing  trends  in  resistivity/ice  saturation 
towards  the  outer  extents  of  the  lake  are  captured  in  the  pdfs,  but  are  subtle.  (5)  Within  the  silt 
layer  at  early  times  before  the  talik  is  fully  through-going  (Figure  5. 5.2. 5  B,  D),  the  AEM  data 
are  insensitive  to  which  layer  is  present,  hence  the  bi-modal  resistivity  distribution  with  peaks 
associated  with  characteristic  silt  and  gravel  values.  This  ambiguity  disappears  at  later  times 
when  the  low-resistivity  unfrozen  silt  layer  extends  to  the  base  of  the  unfrozen  gravels,  which  is 
a  more  resolvable  target  (Figure  5. 5.2. 5  F). 
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Figure  5,5,2.5  Resistivity  probability  distributions  for  tbe  hydrostatic  6-m-deep  lake 
scenario  at  simulation  times  100  years  (A-B),  240  years  (C-D),  and  1,000  years  (E-F), 
Shading  in  each  image  represents  the  probability  distribution  at  depths  of  15  m  (A,  C,  E) 
and  50  m  (B,  D,  F)  from  the  lake  center  (r  =  0  m)  to  the  edge  of  the  model  (r  =  1800  m). 
Dashed  lines  indicate  the  true  resistivity  values. 
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The  final  component  of  this  analysis  focused  on  the  upper  sand  layer  (Unit  1),  which  is  generally 
too  thin  (2  m)  and  resistive  (>  600  ohm-m)  to  be  resolved  using  AEM  data;  though  may  be 
imaged  using  other  ground-based  electrical  or  electromagnetic  geophysical  methods.  Seasonal 
thaw  and  surface  runoff  causes  locally  reduced  resistivity  values  in  the  upper  1  m,  which  is  still 
too  shallow  to  resolve  adequately  using  AEM  data.  In  practice,  shallow  thaw  and  sporadic 
permafrost  trends  are  observed  to  greater  depths  in  many  locations,  including  inactive  or 
abandoned  channels.  To  simulate  these  types  of  features,  the  shallow  resistivity  structure  of  the 
6  m-deep  hydrostatic  lake  scenario  at  1,000  years  was  manually  modified  to  include  three 
synthetic  ‘channels’.  These  channels  are  not  intended  to  represent  realistic  pathways  relative  to 
the  lake  and  the  hydrologic  simulations;  they  are  solely  for  the  purpose  of  illustrating  the  ability 
to  resolve  shallow  resistivity  features. 

Eigure  5. 5. 2. 6  shows  the  three  channels  in  a  zoomed-in  view  of  the  uppermost  portion  of  the 
model  outside  the  lake  extent.  Each  channel  is  100  m  wide,  but  with  different  depths:  1  m  (half 
the  Unit  1  thickness),  2  m  (full  Unit  1  thickness)  and  3  m  (extending  into  the  top  of  Unit  2). 
Analysis  of  AEM  data  simulated  for  this  model,  presented  as  the  McMC  average  model,  are 
shown  in  Eigure  5. 5. 2. 6  B.  All  three  channels  are  clearly  identified,  but  their  thicknesses  and 
resistivity  values  are  overestimated  and  cannot  be  distinguished  from  one  another.  To  explore 
the  possibility  of  better  resolving  these  shallow  features,  synthetic  EM  data  were  simulated  using 
the  characteristics  of  a  ground-based  multi-frequency  EM  tool  (the  GEM-2  instrument)  that  can 
be  hand  carried  or  towed  behind  a  vehicle,  and  was  commonly  used  the  shallow  investigations  of 
the  larger  project.  The  McMC  average  model  result  for  the  simulated  shallow  EM  data  is  shown 
in  Eigure  5. 5. 2. 6  C.  Channel  thicknesses  and  resistivity  values  are  better  resolved  compared  with 
the  AEM  result,  though  the  1  m-deep  channel  near  r  =  800  m  appears  both  too  thick  and  too 
resistive.  In  addition,  the  shallow  EM  data  show  some  sensitivity  to  the  interface  at  2-m  depth 
between  frozen  silty  sands  and  frozen  gravels,  though  the  depth  of  this  interface  is  over¬ 
estimated  due  to  the  limited  sensitivity  to  these  very  resistive  features. 

In  the  specific  examples  of  lake  talik  evolution  presented  here,  which  were  modeled  after  the 
physical  setting  of  the  Yukon  Elats,  Alaska,  AEM  data  were  shown  to  be  generally  capable  of 
resolving  large-scale  permafrost  and  geological  features,  as  well  as  thermally  and  hydrologically 
induced  changes  in  permafrost  over  time.  The  Bayesian  McMC  analysis  provided  useful  details 
about  model  resolution  and  uncertainty  that  cannot  be  assessed  using  traditional  inversion 
methods  that  produce  a  single  ‘best’  model.  A  fortuitous  aspect  of  the  Yukon  Elats  model  is  the 
fact  that  the  silt  layer  (Unit  3)  is  relatively  conductive  compared  with  the  overlying  gravels  (Unit 
2),  making  it  a  good  target  for  electromagnetic  methods.  If  the  order  of  these  layers  were 
reversed,  or  if  the  base  of  permafrost  were  hosted  in  a  relatively  resistive  lithology,  AEM  data 
would  not  likely  resolve  the  overall  structure  with  such  good  fidelity.  In  addition,  knowledge  of 
the  stratigraphy  helps  to  remove  the  ambiguity  between  unfrozen  gravels  and  frozen  silts,  which 
have  similar  intermediate  resistivity  values.  The  methods  developed  here  that  use  a  physical 
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property  model  to  link  hydrologie  and  geophysieal  properties  provide  the  neeessary  framework 
to  test  other  more  challenging  hydrogeological  scenarios. 
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Figure  5,5,2,6  Comparison  of  airborne  and  ground-based  measurements  for  recovering 
shallow  tbaw  features.  True  shallow  resistivity  structure  extracted  from  the  hydrostatic  6- 
m-deep  lake  scenario  at  a  simulation  time  of  1,000  years,  shown  outside  of  the  lake  extent 
(distance  >  500  m)  (A),  Three  shallow  low-resistivity  channels  with  thicknesses  1  m,  2  m, 
and  3  m  were  added  to  the  resistivity  model  to  provide  added  contrast,  McMC-derived 
results  using  simulated  AEM  data  (B)  and  ground-based  EM  data  (C)  illustrate  the 
capability  of  these  systems  to  image  shallow  features. 

Analysis  of  AEM  surveys  provide  a  means  for  remotely  detecting  subsurface  electrical  resistivity 
associated  with  the  co-evolution  of  permafrost  and  hydrologic  systems  over  areas  relevant  to 
catchment-scale  and  larger  processes.  Coupled  hydrogeophysical  simulations  using  a  physical 
property  relationship  that  accounts  for  the  effects  of  lithology,  ice  saturation,  and  temperature  on 
electrical  resistivity  provide  a  systematic  framework  for  exploring  the  geophysical  response  to 
various  scenarios  of  permafrost  evolution  under  different  hydrological  forcing.  This  modeling 
approach  provides  a  means  of  robustly  testing  the  interpretation  of  AEM  data  given  the  paucity 
of  deep  boreholes  and  other  ground  truth  data  that  are  needed  to  characterize  subsurface 
permafrost.  A  robust  uncertainty  analysis  of  the  geophysical  simulations  provides  quantitative 
information  about  what  types  of  features  can  be  resolved  using  AEM  data  given  the  inherent 
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resolution  limitations  of  geophysical  measurements  and  ambiguities  in  the  physical  property 
relationships.  In  the  scenarios  considered  here,  we  have  shown  that  large-scale  geologic  and 
permafrost  structure  is  accurately  estimated.  Sublacustrine  thaw  can  also  be  identified,  but  the 
specific  geometry  of  partial  ice  saturation  beneath  lakes  can  be  poorly  resolved  by  AEM  data. 
Understanding  the  geophysical  response  to  known  simulations  is  helpful  both  for  guiding  the 
interpretation  of  existing  AEM  data,  and  also  to  plan  future  surveys  and  other  ground-based  data 
acquisition  efforts. 


5,6  Permafrost  Response  to  Lake  Recession 

Research  at  Twelvemile  Eake  in  the  Yukon  Elats  of  Alaska  revealed  previously  unknown  details 
about  permafrost  distribution  and  dynamics  associated  with  lake  recession  in  discontinuous 
permafrost.  Airborne  electromagnetic  surveys  suggest  that  unfrozen  conditions  exist  below 
Twelvemile  Eake  (i.e.,  an  open  talik)  and  relatively  thick  (--90  m)  permafrost  exists  outside  the 
perimeter  of  the  highstand  of  the  lake  in  the  mid-1980s.  Therefore,  as  the  lake  receded,  exposed 
open  meadow  was,  at  least  initially,  free  of  permafrost.  Over  the  past  two  decades  willow  shrubs 
have  grown  in  discrete  clumps/bands  throughout  the  previously  open  meadow  area  through  the 
process  of  ecological  succession.  Ground-based  electrical  geophysical  surveys  conducted  as  part 
of  this  project  in  late  summer  2011  and  2012  along  a  transect  of  mixed  open-meadow  and  willow 
shrub  showed  shallow  zones  with  permafrost  located  proximal  to  maturing  shrubs,  while  more 
open  areas  did  not  show  evidence  of  frozen  ground.  Manual  frost  and  temperature  probing 
confirmed  these  results.  A  modeling  approach  was  applied  to  provide  mechanistic  understanding 
into  the  observations  of  permafrost  aggradation  in  Twelvemile  Eake  receding  margin. 

5.6.1  Materials  and  Methods 

A  suite  of  100-year  vertical  one-dimensional  (1-D)  computer  simulations  was  used  to  represent 
various  ecosystem  effects  (reduced  recharge  and  summer  cooling  from  shading)  and  climate 
change  using  the  modified  SUTRA  code  to  simulate  freeze-thaw.  Unsaturated  zone  dynamics 
had  not  been  considered  in  previous  numerical  modeling  of  permafrost  aggradation.  The 
simulations  were  designed  to  (1)  test  the  hypothesis  that  permafrost  can  form  under  current 
conditions  in  response  to  lake  recession,  (2)  gain  insight  into  factors  controlling  permafrost 
formation,  and,  (3)  evaluate  the  longer-time  response  to  climate  warming.  The  open  meadow 
with  summer  recharge  (OM-R)  scenario  serves  as  a  base  case  representing  receded  lake 
conditions  prior  to  shrub  succession.  Prescribed  model  ground-surface  temperature  and  recharge 
boundary  conditions  for  OM-R  reflect  local  observations  and  history  (NRCS  SNOw  TEEemetry 
Station  #961).  Eor  comparison,  we  considered  various  combinations  of  shading  at  two  levels 
(Shi  and  Sh2),  summer  recharge  at  the  current  level  and  zero  level  (R  and  OR,  respectively),  and 
climate  warming  (CW),  for  a  total  of  7  models:  OM-R,  Shl-R,  Sh2-R,  OM-OR,  Shl-OR,  Sh2-0R, 
and  Sh2-0R-CW.  Eor  the  shading  scenarios,  moderate  reductions  of  1  °C  (Shi)  and  2  °C  (Sh2)  in 
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peak  summer  shallow  soil  temperatures  (defined  at  0.06  m  model  depth)  were  based  on  observed 
willow  shrub  shading  effeets  in  a  similar  setting  (Smith,  1975);  these  values  approximate  the 
average  apparent  eooling  effeet  of  0.6  °C  at  0.5  m  depth  under  diserete  willow  elumps  (n=5) 
eompared  to  adjaeent  open  meadow  areas  (n=5)  measured  at  the  field  site  in  late  summer  2013. 

Summer  reeharge  averages  0.1  m  in  the  Twelvemile  Lake  area  (SNOTEL  #961),  well  below  the 
maximum  potential  for  willow  evapotranspiration;  therefore,  effeetive  zero-recharge  (OR)  is  a 
plausible  condition.  The  strongest  ecosystem  effect  considered  was  a  2  °C  reduction  in 
maximum  shallow  soil  temperatures  and  no  summer  recharge  (Sh2-0R).  The  silt-loam  soils  at 
Twelvemile  Lake  were  simulated  using  the  recently  measured  unsaturated  soil-water  freezing 
and  soil-water  retention  characteristics  of  similar  materials  (Watanabe  et  ah,  2011;  Kurylyk  and 
Watanabe,  2013). 

The  domain  spanned  an  elevation  of  +1  m  to  a  depth  of  20  m,  with  the  datum  being  the  ground 
surface.  The  upper  1  m  represented  a  thermal  boundary  layer  with  a  sinusoidal  annual  thermal 
input.  The  upper  model  domain  spanned  0-5  m  in  depth  and  consisted  of  0.0 1-m  vertical 
elements.  The  water  table  at  -3.0  m  was  represented  as  a  zero-pressure  boundary,  simulating  the 
control  of  a  constant  lake  level  on  water  table,  based  on  observations  over  the  2011-2012  field 
seasons.  Therefore,  the  model  domain  from  0-3  m  depth  was  variably  saturated  through  time. 
The  lower  model  domain,  from  5-20  m  depth,  comprised  a  vertical  column  of  elements  0.1-m 
thick  to  simulate  the  deeper  saturated  zone  and  the  geothermal  gradient  from  below  the 
penetration  depth  of  the  annual  temperature  envelope. 

The  initial  condition  for  each  model  was  the  result  of  an  80-yr  equilibration  run  of  the  OM-R 
model  using  a  sinusoidal  annual  thermal  input  with  amplitude  of  16  °C  and  a  mean  of  0.5  °C, 
similar  to  observed  ground  temperatures  within  the  dried  lake  boundary  near  the  study  site.  The 
phase  of  the  annual  temperature  sinusoid  was  shifted  forward  in  time  by  3  months  (7i/2)  to  best 
coincide  the  onset  of  snowmelt  with  a  few-day  lag  of  boundary  temperature  rising  above  0  °C  in 
April.  Once  initial  conditions  were  established,  predictive  models  were  run  for  100-yr 
simulations,  incorporating  varied  effects  of  willow  shading  by  reducing  annual  temperature 
amplitude  and  mean  by  similar  magnitudes;  in  this  way  the  summer  thermal  maximum  is 
reduced  but  the  cold  winter  minimum  is  preserved.  Note  that  potential  winter  insolation  effects 
of  trapped  snow  were  not  incorporated. 

The  climate  warming  scenario  was  created  by  applying  a  trend  to  the  thermal  boundary  that 
resulted  in  a  +2.55  °C  change  in  mean  temperature  at  0.06  m  depth  after  100  yr;  the  2.55  °C 
value  was  chosen  by  applying  a  conversion  factor  of  0.85  for  “low  shrub-scrub”  land-cover 
(Jorgenson  and  Kreig,  1988)  on  soils  to  the  moderate  Scenarios  Network  for  Alaska  and  Arctic 
Planning  (SNAP)  (Scenarios  Network  for  Alaska  and  Arctic  Planning,  University  of  Alaska, 
2013)  predicted  air  warming  over  the  simulated  time  period.  Little  regional  trend  in  recharge  is 
predicted  by  SNAP,  so  the  recharge  boundary  was  left  unchanged.  Once  the  +2.55  °C 
background  temperature  trend  was  established,  the  climate  warming  model  scenario  was  re-run 
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with  the  added  “strong  willow  effeet”  of  2  °C  summer  shading  and  no  summer  reeharge  (model 
Sh2-0R-CW). 

A  time  varying  reeharge  boundary  was  set  at  the  ground  surfaee  to  represent  three  distinet 
phases:  (1)  Winter,  no  recharge  (ordinal  days  0-100  and  275-365),  (2)  Spring  snowmelt  recharge 
(days  100-120),  (3)  Summer  recharge  (days  120-275).  Meltwater  temperature  was  0.5  °C,  and 
volume  was  determined  from  average  Fort  Yukon  data  (SNOTEL  #961)  as  0.47  m  snow  with  0.2 
snow/water  equivalent.  The  summer  recharge  phase  delivers  0.1  m  of  precipitation;  this 
recharge  was  assigned  a  value  of  12  °C  to  represent  observed  area  conditions,  although  models 
were  found  to  be  relatively  insensitive  to  this  assumed  temperature.  To  simulate  the  effects  of 
transpiration/interception  by  the  willow  shrub,  summer  recharge  was  set  to  0  m  (OR  models), 
while  the  melt  pulse  was  preserved  as  it  was  assumed  the  shrubs  are  not  foliated  and  therefore 
not  using  soil  water  at  that  time. 

Various  unsaturated  soil  properties  were  investigated,  and  simulations  were  found  to  be  strongly 
sensitive  to  both  the  freezing  function  and  residual  saturation  at  strong  negative  (drying) 
pressure.  Therefore,  great  care  was  taken  to  match  model  parameters  with  the  general  soil  type 
observed  in  the  field.  A  single  homogeneous  soil  was  used  for  all  models  and  based  on  a 
classification  of  “silt-loam”  and  observations  from  soil  pits.  The  unsaturated  soil-water  freezing 
and  soil-water  retention  characteristics  for  this  soil  type  were  based  on  a  recent  laboratory 
investigation  (Watanabe  et  ah,  2011),  as  was  the  silt-loam  porosity  of  0.46.  The  soil-water 
freeze  function  (relating  temperature  to  ice  content)  was  simulated  with  a  4-section  piecewise- 
linear  function  with  a  minimum  residual  unfrozen  water  content  of  16.6%  (of  total  porosity)  at  - 
20  °C.  The  soil-water  retention  curve  was  simulated  with  a  piecewise-  linear  function  with 
residual  saturation  at  50,000  kPa  of  8%.  The  soil  had  a  base  permeability  of  1x10'  m  ,  with  a 
frozen  minimum  permeability  of  1x10'  m  .  The  relatively  small  reduction  in  permeability 
when  frozen  was  based  on  observation  of  water  flow  through  frozen  silt-loam,  compared  to  other 
frozen  soils  that  do  not  permit  flow  when  frozen  (Watanabe  et  ah,  2011,  2013);  permeabilities 
less  than  approximately  1x10'  m  essentially  function  as  confining  material  within  the  model 
and  could  have  induced  surface  ponding,  a  condition  for  which  there  was  no  field  evidence. 

5.6.2  Results  and  Discussion 

Model  scenarios  presented  and  discussed  in  this  section  include: 

(1)  OM-R:  Open  meadow  base  case  model  with  seasonal  recharge  as  specified  in  section 
5.6.1 

(2)  OM-OR:  Open  meadow  model  with  no  seasonal  recharge 

(3)  Shl-R:  Moderate  shading  effect  model  with  1  °C  summer  cooling  and  seasonal  recharge 

(4)  Shi -OR:  Moderate  shading  effect  model  with  1  °C  summer  cooling  and  no  seasonal 
recharge 

(5)  Sh2-R:  Strong  shading  effect  model  with  2  °C  summer  cooling  and  seasonal  recharge 
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(6)  Sh2-0R;  Strong  shading  effect  model  with  2  °C  summer  cooling  and  no  seasonal  recharge 

(7)  Sh2-0R-CW;  Strong  shading  effect  model  with  2  °C  summer  cooling,  no  seasonal 
recharge,  and  linear  climate  warming  trend  superimposed  on  seasonal  temperature  cycle 

Results  from  the  OM-R  base  model  show  seasonal  frost  and  temperature  propagation  through 
time  (Figure  5.6.2. 1).  Annual  frost  and  temperature  dynamics  define  envelopes  that  encompass 
the  range  of  simulated  conditions  with  depth  along  the  profile;  for  OM-R,  the  maximum  and 
mean  shallow  soil  temperatures  are  12.3  and  0.6  °C,  respectively.  Under  the  expected  ambient 
conditions  of  the  open  meadow  dried  lakebed,  appreciable  permafrost  does  not  form.  Notably, 
annual  downward  propagation  of  ice  from  the  previous  winter  through  the  soil  peaks  near  the 
water  table  in  early  September,  yet  ice  content  drops  to  zero  for  the  rest  of  the  annual  cycle 
(Figure  5.6.2. 1  a).  This  strong  phase  shift  has  important  implications  for  field  interpretations  of 
frozen  ground.  The  presence  of  low  ice-content  deep  frost  late  in  the  summer  may  indicate 
seasonally  frozen  ground  rather  than  permafrost.  The  simulated  low  year-round  ice  saturation 
(3%  of  porosity)  at  the  water  table  in  OM-R  (Figures  5.6.2.1a,  5.6.2.2a)  agrees  with  site 
observations  from  a  soil  pit  dug  in  open  meadow,  where  stiff  soil  and  temperatures  close  to  0  °C 
were  detected  in  a  discrete  band  just  above  the  water  table,  which  supports  the  OM-R  model. 
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Figure  5,6,2, 1  Three-year  ice  saturation  time  series  for  multiple  depths  in  the  open-meadow 
model  showing  that  ice  saturation  peaks  around  1-m  depth,  and  propagates  downward 
while  shifting  forward  in  phase  with  yearly  peak  ice  content  indicated  hy  red  dots;  the 
starred  dot  corresponds  to  the  time  year  of  2012  data  collection  (a).  Smoothed  (1  year 
average)  time  series  of  ice  saturation  at  2,5  m  depth  in  all  models,  showing  the  accelerated 
frost  formation  caused  hy  uptake  of  summer  recharge  and  permafrost  aggradation  and 
subsequent  degradation  under  climate  change. 

Simulation  OM-OR  indicates  that  willow  shrub  uptake  of  summer  recharge  does  not  result  in 
appreciable  permafrost  aggradation  within  the  Twelvemile  Lake  basin.  When  recharge  is 
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removed  from  the  baseline  model,  the  original  diserete  band  of  lee-saturation  at  the  water  table 
interface  is  slightly  enhanced  to  7%;  hence  shrub-driven  reduction  in  recharge  cannot  explain  the 
thick,  shallow  sequence  of  frost  observed  in  the  field  beneath  shrub  areas.  Recharge  reduction 
has  a  much  stronger  effect  when  combined  with  summer  shading.  Figure  5.6.2.1b  shows  how 
frost  development  at  greater  depth  in  the  unsaturated  zone,  characterized  by  multiple  model  time- 
series  at  2.5  m  depth,  is  greatly  accelerated  in  the  Shi  model  when  summer  recharge  is  removed 
in  Shi -OR.  There  is  muted  enhancement  of  freezing  in  the  stronger  shading  Sh2-0R  model, 
compared  to  Sh2-R,  indicating  the  importance  of  unsaturated  zone  heat  advection  to  permafrost 
aggradation  peaks  under  moderate  shading  conditions.  Thus,  the  uptake  of  summer  recharge  and 
subsequent  reduction  of  infiltration  alone  does  not  result  in  strong  permafrost  aggradation,  but 
when  combined  with  summer  shading,  recharge  reduction  does  result  in  appreciable 
enhancement  of  freezing,  particularly  in  the  Shi -OR  model. 

Thick  sequences  of  new  shallow  permafrost  accumulate  in  models  Shi -OR  (5.7  m  thick)  and 
Sh2-0R  (7.7  m  thick),  which  have  both  shading  and  reduced  recharge,  as  shown  in  Figure  5. 6.2. 2 
b,c,  respectively.  Dynamic  steady-  state  permafrost  is  0.2  m  shallower  in  model  Sh2-0R  at  1.6  m 
(Figure  5.6.2.2c  and  5.6.2.3a);  this  depth  is  within  0.16  m  of  that  observed  near  willow  shrub  in 
the  field.  Much  of  the  frost  accumulation  takes  place  in  the  unsaturated  zone  in  the  first  25  yr, 
particularly  with  increased  shading  (Sh2-0R).  As  the  unsaturated  zone  permafrost  quasi- 
equilibrates,  substantial  freezing  of  the  saturated  zone  occurs  over  the  remainder  of  the 
simulations.  The  Shi -OR  and  Sh2-0R  models  are  most  dissimilar  in  quasi-equilibration  time  and 
saturated  zone  freezing  dynamics.  Ice  saturation  in  the  shallow  aquifer  of  both  models  exceeds 
50%  of  the  available  porosity,  thus  restricting  flow  and  the  advection  of  heat.  This  result  shows 
how  critical  yet  subtle  differences  in  shrub  growth  and  subsequent  shading  may  have  pronounced 
impacts  on  shallow  groundwater  flow  and  the  resultant  lake  water  budget. 
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Figure  5,6,2,2  Ice  saturation  envelope  encompassing  annual  frost  variation  (including  3 
date-specific  profiles)  and  aggraded  permafrost  ice  saturation  after  99  yr  of  simulation  for 
models:  open-meadow  (OM-R)  (a),  moderate  shrub  effect  (Shl-OR)  (b),  strong  shrub  effect 
(Sh2-0R)  (c),  and  strong  shrub  effect  with  climate  change  (Sh2-0R-CW)  (d). 

The  SNAP-guided  climate  warming  (CW)  trend  was  imposed  on  the  hydrologic  scenario  that 
yielded  the  thickest  aggraded  permafrost  under  current  conditions  (Sh2-0R).  The  Sh2-0R-CW 
simulation  (Figure  5.6.2.3b)  shows  initial  permafrost  aggradation  on  the  scale  of  several  meters 
in  thickness,  similar  to  Sh2-0R,  though  after  approximately  20  yr  the  simulations  diverge  and 
permafrost  begins  to  degrade  for  Sh2-0R-CW.  Ice  saturation  peaks  at  45%  with  climate 
warming  and  60%  without  warming.  The  maximum  depth  reached  by  permafrost  subjected  to 
climate  warming  is  6.1  m  (9.3  m  with  no  warming),  and  this  maximum  occurs  at  year  45,  after 
which  permafrost  degrades.  Under  climate  warming,  frozen  ground  first  thaws  downwards  in  the 
unsaturated  zone,  vanishing  by  year  69  (Figure  5.6.2.1b).  Along  with  permafrost  degradation, 
the  seasonally  active  freeze/thaw  layer  thins  from  a  zone  extending  from  land  surface  to  below 
the  water  table  to  a  zone  only  reaching  1.8  m  depth  by  year  100  (Figure  5.6.2.2d).  This  thinning 
also  is  accompanied  by  a  broadening  in  the  annual  temperature  envelope  to  incorporate 
temperatures  at  depths  that  are  warmer  than  observed  in  any  other  model);  the  envelope 
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converges  on  the  geothermal  gradient  at  several  degrees  above  zero  approximately  3  m  below  the 
water  table.  The  effects  of  expected  air  temperature  inereases  on  near-surface  soil  temperatures 
may  be  countered,  at  least  to  some  degree,  by  reductions  in  the  thickness  and  duration  of 
snowpaek  and  partieularly  by  the  later  timing  of  initial  snow  aceumulation.  The  modeling 
approach  presented  here  to  represent  elimate  ehange  assumes  such  effects  of  potential  changes  in 
snowpaek  are  negligible  or  at  least  overshadowed  by  the  influence  of  warming  air  temperatures 
in  the  study  region. 


time  (yr)  time  (yr) 


Figure  5,6,2,3  Strong  shrub  effect  model  ice  saturation  (Sh2-0R)  (a),  and  shrub  effect 
reduced  by  climate  change  model  ice  saturation  (Sh2-0R-CW),  Dashed  line  indicates  static 
water  table.  Ice  saturation  is  defined  as  fraction  of  pore  space  filled  with  frozen  water. 

Through  numerical  simulation,  this  study  revealed  a  mechanism  to  explain  observed  permafrost 
aggradation  under  eurrent  conditions.  Permafrost  aggradation  following  lake  recession  and 
vegetation  succession  may  be  an  important  eomponent  in  the  eycle  of  lake  or  wetland  contraetion 
and  expansion.  Permafrost  formation  alters  the  surfaee-water/groundwater  exehange  capacity 
and  near-surface  fluxes  of  water  and  solutes.  Permafrost  aggradation  may  act  to  reverse  the  net 
negative  water  flux  in  a  reeeding  lake,  thereby  promoting  the  natural  eycle  of  lake-level 
fluetuations  observed  in  lake  eores  in  the  region.  Should  warming  oceur  as  predieted,  however, 
simulation  results  point  to  termination  of  the  lake-retreat/ecosystem  mechanism  as  a  driver  of 
permafrost  aggradation  before  the  turn  of  the  next  century. 
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5,7  Regional  Synthesis  of  Yukon  Flats  Lake  Dynamics  and  Linkages  to  Permafrost 
Conditions 

The  multi-scale  analyses  of  permafrost  and  hydrology  interplay  together  with  the  multi-scale 
characterization  of  permafrost  in  the  lake-rich  landscape  of  the  Yukon  Flats  provided  a 

foundation  for  examining  linkages  between  permafrost  distribution  and  lake  surface-area  change 

2 

in  cold  regions.  Fiere,  a  large-scale  examination  of  these  linkages  was  made  over  a  5,150-km 
area  of  Yukon  Flats  by  evaluating  the  relationship  between  lake  surface-area  changes  during 
1979-2009,  derived  from  Landsat  satellite  data  and  sublacustrine  groundwater- flowpath 
connectivity  inferred  from  the  AEM  survey  of  permafrost  (Figure  5.7.0. 1).  Furthermore,  insight 
gained  from  our  modified  SUTRA  modeling  analyses  was  used  to  elucidate  empirical  results 
generated  by  the  data  comparison  analysis. 
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^  Expanding  lake  for  surface  area  change 


Figure  5,7,0, 1  Location  of  study  area  (a,  c),  and  lake  surface-area  change  hy  physiographic 
unit  (h).  The  upper  and  lower  numbers  in  the  har  charts  of  (h)  are  the  percentages  of  total 
lakes  (~8500)  and  total  lake  surface  area  (53,6  x  10^  ha)  in  the  study  area. 
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5.7.1  Materials  and  Methods 


Lake  taliks  (totaling  153)  were  eharaeterized  from  eleetrieal  resistivity  eross-seetions  along  392 
flight-line  km  of  the  AEM  survey  (Figure  5.7.0. 1).  Interpretation  of  the  resistivity  eross-seetions 
followed  that  of  Minsley  et  al.  (2012).  The  resistivity  of  geologie  materials  varies  with  both 
roek/soil  type  and  thaw  state  (i.e.,  frozen  versus  unfrozen)  (Hoekstra  et  ah,  1975),  potentially 
resulting  in  an  overlap  in  resistivity  between  different  materials  having  eontrasting  thaw  states. 
To  reduee  this  ambiguity,  we  ineorporated  knowledge  about  the  known  depositional  history  of 
the  Yukon  Flats  (Clark  et  al.,  2009;  Williams,  1962)  to  develop  a  lithologie  model  and  place 
constraints  on  the  spatial  distribution  of  different  materials.  The  lithologic  model  consists  of  2 
layers:  an  upper  layer  of  fluvial  gravel  ('-'15-50  m  thick)  and  an  underlying  layer  of  lacustrine 
silt.  The  observed  resistivity  transition  occurring  at  the  interface  between  frozen  gravel  and 
frozen  silt,  and  between  frozen  silt  and  unfrozen  silt,  in  a  deep  borehole  in  Fort  Yukon  (Clark  et 
ah,  2009)  (transitions  illustrated  along  line  F-F’  of  Figure  5.7.1.1a)  were  used  to  constrain  our 
model  relating  thaw  state  to  resistivity  for  each  lithology  (Figure  5. 7. 1.1,  bottom).  The  interface 
between  fluvial  gravel  and  underlying  lacustrine  silt,  given  their  depositional  environment,  is 
expected  to  be  horizontal  at  the  scale  of  individual  lakes.  Therefore,  at  a  given  depth,  it  is 
reasonable  to  assume  that  lateral  resistivity  transitions  across  lake  basins  are  related  to  changes  in 
thaw  state,  rather  than  material  type.  Shallow  groundtruth  data  and  a  discussion  about  possible 
thaw  state  interpretation  errors  are  provided  as  auxiliary  material. 

Sublacustrine  silt  is  considered  to  be  unfrozen  if  the  vertical  distance  between  the  gravel-silt 
contact  and  interpreted  bottom  of  permafrost,  tf,  is  less  than  a  threshold  value,  set  to  10  m  to 
account  for  the  uncertainty  in  interface  positions  (Figure  5.7. 1.1).  One  of  four  possible 
combinations  of  thaw  states  for  sublacustrine  gravel  and  silt  (Type  A-D  cases.  Figure  5. 7. 1.1)  are 
assigned  to  each  lake.  Type  A  cases  are  assumed  to  represent  open  taliks,  and  Type  D  cases  are 
assumed  to  represent  closed  taliks  with  a  phase  boundary  occurring  near  the  gravel-silt  contact. 
In  Type  B  and  C  cases,  the  gravel  is  assumed  to  be  completely  or  partially  frozen. 

Determination  of  shrinking,  stable,  and  expanding  lakes  was  based  on  the  regression  analysis  by 
Rover  et  al.  (2012)  using  Fandsat  satellite  data  (USGS,  2011).  This  analysis  used  decision  tree 
classification  models  (average  model  accuracy  98.9%)  to  map  lake  surface  areas,  and  linear 
regression  to  determine  if  lake  surface  areas  have  decreasing  (“shrinking”  lakes),  no  (“stable” 
lakes),  or  increasing  (“expanding”  lakes)  trends  at  the  95%  level  of  statistical  significance  (p- 
value  <  0.05)  for  the  period  of  1979-2009.  These  surface  area  trends  were  calculated  from 
available  cloud  free  observations  during  May-September  (17-20  observations  per  lake). 
Existence  of  a  statistical  association  between  lake  surface-area  trend  and  sublacustrine  thaw  state 
was  determined  using  contingency  tables  with  chi-square  tests  of  independence  (provided  as 
auxiliary  material).  Fake  surface-area  trend  was  considered  to  be  “significantly  associated”  with 
sublacustrine  thaw  state  if  the  null  hypothesis,  postulating  that  area  trend  and  thaw  state  are 
independent,  was  rejected  at  the  95%  level  of  significance  (p-value  <  0.05). 
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Figure  5,7,1, 1  Examples  of  the  four  different  gravel-silt  thaw  states  (Types  A-D) 
interpreted  from  electrical  resistivity  cross  sections  of  the  airborne  electromagnetic  survey. 
The  vertical  distance  between  the  gravel-silt  contact  and  the  interpreted  bottom  of 
permafrost  tf,  is  used  to  determine  whether  silt  if  unfrozen  (tf  <  10  m)  or  frozen  (tf  >10  m). 
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5.7.2  Results  and  Discussion 


The  occurrence  of  open  taliks  (Type  A)  is  not  significantly  associated  with  lake  area  trend  (p- 
value  =  0.37),  as  illustrated  by  the  similar  percentages  of  area  trends  for  Type  A  cases  and  all 
lakes  (Figure  5.7.2.1a).  This  suggests  that  the  formation  of  open  taliks  is  not  a  primary 
mechanism  of  the  observed  lake  surface-area  dynamics  in  the  Yukon  Flats.  The  lack  of 
association  between  open  taliks  and  lake  area  trends  also  reduces  the  likelihood  that  changes  in 
regional  groundwater  flow,  which  influence  lakewater  /  deep-groundwater  exchange  through 
sublacustrine  open-taliks,  account  for  the  observed  lake  dynamics. 

In  contrast  to  open  taliks,  the  occurrence  of  unfrozen  sublacustrine  gravel  is  significantly 
associated  with  lake  area  trend  (p-value  =  0.04).  For  example,  lakes  overlying  unfrozen  gravel 
are  2.5  times  as  likely  to  be  shrinking  as  lakes  not  overlying  unfrozen  gravel  (Figure  5.7.2.1b). 
Lakes  overlying  unfrozen  gravel  and  frozen  silt  (Type  D  cases),  most  prevalent  in  the  Loess  Base 
(Figure  5.7.2.1c),  are  particularly  susceptible  to  shrinkage  (Figure  5.7.2.1a).  These  observations 
suggest  that  changes  in  lakewater  /  groundwater  exchange  as  a  potential  driver  of  lake  volume 
evolution  result  more  likely  from  shallow  (few  tens  of  meters),  rather  than  deeper  (-'50-100  m), 
thermal  changes  in  permafrost.  Shallow  thermal  changes  influencing  lakewater  /  groundwater 
exchange  would  likely  need  to  occur  in  terrestrial  areas  of  watersheds  in  order  to  provide  lateral 
groundwater-fiowpath  connectivity  to  and  from  sublacustrine  aquifers,  possibly  including 
deepening  of  the  permafrost  table  and  growth  of  supra-  and/or  intra-permafrost  taliks.  Such 
processes  could  allow  substantial  groundwater  exchange  between  watersheds  in  the  study  area 
owing  to  the  high  hydraulic  conductivity  of  unfrozen  gravel.  In  addition,  the  occurrence  of 
unfrozen  gravel  is  more  spatially  variable  across  the  landscape  than  deeper  silt,  as  indicated  by 
the  greater  standard  deviation  of  unfrozen  gravel  than  deeper  unfrozen  silt  across  physiographic 
units  (19%  in  Figure  5. 7. 2. Id  versus  13%  in  Figure  5.7.2. le). 

The  occurrence  of  frozen  sublacustrine  silt  is  significantly  associated  with  lakes  that  are 
expanding  (p-value  =  0.05).  Lakes  overlying  frozen  silt  are  2.7  times  as  likely  to  be  expanding  as 
lakes  overlying  unfrozen  silt  (Figure  5. 7. 2. If).  Given  the  low  hydraulic  conductivity  of  frozen 
silt,  upwelling  of  deep  groundwater  is  not  likely  to  be  a  significant  source  of  recharge  to  many  of 
these  expanding  lakes.  Rather,  enhanced  water  supply  to  these  lakes  may  follow  shallow 
flowpaths  through  unfrozen  gravel,  or  overland  fiowpaths  from  adjacent  water  bodies.  The 
presence  of  frozen  silt  below  many  of  the  expanding  lakes  may  be  an  indication  regarding  their 
age  and/or  persistence  of  being  filled  with  water  -  many  of  them  may  be  too  young  to  have 
formed  open  taliks,  or  may  undergo  cycles  of  filling  and  drainage  thereby  allowing  their  lake 
bottoms  to  freeze  intermittently.  The  observed  high  occurrence  of  expanding  lakes  along  rivers 
and  creeks  (Figure  5.7.0.1b),  coupled  with  the  possibility  that  many  of  these  lakes  are  young, 
could  be  an  indication  of  recent  increases  in  regional  groundwater  discharge  through  taliks  in 
low-lying  river  corridors  (Walvoord  et  ah,  2012). 
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(a)  Surface  area  trend: 


(b)  Surface  area  trend: 


(c)  Gravel-silt  thaw  state: 


□  No  change  □  No  change  □Type  A  I  I  Type  B 

□  Shrinking  H  Expanding  □Shrinking  ■  Expanding  ■TypeC  ^ Type  D 


(d)  Gravel  thaw  state: 

I  I  Uncertain 
I  I  Unfrozen  ■  Frozen 


Physiographic  unit 


(e)  Siit  thaw  state: 


Physiographic  unit 


(f)  Surface  area  trend: 

I  I  No  change 
I  I  Shrinking  ■  Expanding 


Silt  thaw  state 


Figure  5,7.2. 1  Relationships  between  lake  surface-area  trends  and  sublacustrine  thaw 
states  (a,  b,  and  f)  and  sublacustrine  thaw  states  and  location  (c-e).  Thaw  states  in  a  and  c 
are  defined  in  Figure  5,7.1. 1,  Numbers  in  bars  of  a,  b,  and  f  are  the  number  of  lakes 
sampled  by  the  AEM  survey.  In  bars  of  the  remaining  panels,  the  upper  and  lower 
numbers  are  the  percentages  of  total  lake  surface  area  from  Rover  et  al.  (2012)  and  the 
number  of  lakes  samples  by  the  AEM  survey,  respectively. 


Results  suggest  that  if  inereased  groundwater  exchange  is  driving  the  observed  lake  shrinkage 
and  expansion,  the  process  likely  involves  shallow  thermal  changes  in  permafrost  (upper  few 
tens  of  meters),  such  as  deepening  of  the  permafrost  table  and  growth  of  supra-  and  intra- 
permafrost  taliks.  In  the  region  studied,  these  key  shallow  aquifers  exhibit  the  highest  hydraulic 
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conductivity  and  greatest  spatial  variability  in  unfrozen  state.  The  thermal  and  hydrologieal 
properties  of  these  aquifers  will  be  partieularly  sensitive  to  elimate  ehange  owing  to  their  elose 
proximity  to  the  atmosphere.  These  regional-seale  results  are  eonsistent  with  previous  studies 
indieating  the  substantial  impaet  of  shallow  permafrost  thaw  on  lakewater-groundwater  exehange 
and  henee  lake  water  budgets.  Reported  thaw  depths  in  these  studies  have  been  up  to  a  few  tens 
of  meters,  generally  beeoming  shallower  northward  where  permafrost  beeomes  eolder  and  more 
eontinuous.  The  integration  of  airborne  geophysies  in  this  study  has  allowed  a  deeper  and  larger- 
seale  inspeetion  of  linkages  between  groundwater-flowpath  eonneetivity  and  ehanges  in  lake 
water  budgets  than  was  available  in  previous  ground-based  studies  alone. 
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6.  Conclusions  and  Implications  for  Future  Research/Implementation 

The  three  main  objeetives  of  this  multi-disciplinary  study  were  successfully  met.  First,  the  USGS 
SUTRA  code  was  modified  to  include  new  general  capabilities  required  to  simulate  the 
hydrologic  consequences  of  permafrost  degradation  and  the  resulting  impact  to  military  and 
federal  infrastructure  and  operations.  Second,  geophysical  investigations  in  an  area  of 
discontinuous  permafrost  of  interior  Alaska  were  performed  to  develop  a  multi-method  approach 
for  mapping  ground  ice  to  be  used  in  any  cold-region  area  and  to  provide  baseline 
characterization  of  the  study  area.  And  lastly,  a  series  of  integrated  numerical  modeling  and 
hydrologic  analyses  informed  by  results  of  the  geophysical  investigations  elucidated  the 
interaction  between  groundwater  and  permafrost  at  scales  ranging  from  site  to  regional. 
Comparison  between  known  cryohydrologic  conditions  predicted  by  hydrologic  modeling  and 
those  inferred  by  geophysical  forward  modeling  provided  a  means  for  supporting  geophysical 
interpretations  and  identifying  remaining  measurement  limitations.  Collectively,  the  tools, 
methods,  and  insight  developed  here  may  be  transferred  to  characterization  and  climate  change 
assessments  of  a  broad  range  of  cold-region  environments  and  applications.  Demonstration  of 
specific  applications  to  the  Yukon  Flats  study  area  is  of  particular  relevance  to  DoD,  since  most 
of  interior  Alaska  DoD  infrastructure  is  located  in  similar  low-lying  regions. 


6,1.  Synthesis  of  Results  and  Research 

Geophysical  measurements  consisting  of  electrical,  seismic,  radar  and  electromagnetic  methods 
all  have  potential  application  for  characterizing  the  distribution  of  permafrost  and  ground  ice, 
with  the  efficacy  and  practicality  of  each  depending  on  site  conditions,  the  depth  of  the  target 
interface  (i.e.,  contact  between  frozen  and  unfrozen  soil),  and  logistical  issues  associated  with 
survey  layout.  The  remoteness  of  some  cold  regions  implies  limited  accessibility  for  carrying  out 
standard  hydrologic  measurements.  For  example,  much  of  the  Yukon  Flats  Basin  of  Alaska  is 
inaccessible  by  roads,  and  transport  of  equipment  must  occur  via  either  river  barge  or  airlift.  In 
this  area  of  Alaska,  installations  of  typical  hydrologic-studies  boreholes  that  penetrate  the  thick 
permafrost  layer  are  expensive  and  present  technical  challenges.  Such  difficulty  in  applying 
standard  approaches  provides  strong  motivation  for  using  geophysical  methods  to  characterize 
the  subsurface  in  such  regions,  because  such  surveys  can  be  accomplished  at  much  lower  cost. 
Airborne  electromagnetic  data  are  most  useful  for  baseline  characterization  of  subsurface 
properties  and  permafrost  distribution  over  large  areas.  Repeat  AEM  surveys  (e.g.  once  every  10 
years)  can  be  valuable  for  monitoring  overall  trends  in  permafrost  distributions  over  large  areas 
as  they  change  resulting  from  varying  climate  and  human  impacts.  Ground-based  geophysical 
surveys  provide  high-resolution  data  sets  that  are  highly  amenable  to  detecting  change  via  repeat 
measurement.  Ground-based  geophysical  methods  are  thus  useful  for  monitoring  subtle  localized 
changes  in  permafrost  near  specific  targets  of  interest  (such  as  roads,  streams,  lakes  and 
buildings).  Regular  ground-based  geophysical  monitoring  is  valuable  for  informing  engineering 
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interventions  for  problems  caused  by  thawing  ground  ice  in  local  target  areas.  Together, 
geophysical  modeling,  thermophysical  hydrologic  modeling,  and  field  observations,  as 
developed  within  this  project,  create  a  synergy  that  provides  greater  insight  into  permafrost- 
hydrology  relations  than  any  individual  approach,  and  can  be  useful  for  future  characterization  of 
coupled  permafrost  and  hydrologic  processes. 

Integrated  hydrogeologic  modeling  and  analysis  to  study  permafrost/groundwater  interaction  has 
been  demonstrated  here  in  a  series  of  studies.  A  groundwater  flow  modeling  analysis  of  the 
Yukon  Flats  Basin  elucidates  the  influence  of  permafrost  on  regional-scale  flow  magnitude  and 
pattern.  At  this  large  scale,  permafrost  degradation  enhances  overall  groundwater  flow,  and  the 
component  of  water  generated  for  streamflow  derived  from  groundwater  increases  relative  to 
runoff  and  shallow  subsurface  pathways.  Presently,  the  Yukon  Flats  Basin  is  coarsely  mapped  as 
spanning  the  continuous-discontinuous  permafrost  transition  that  model  analysis  shows  to  be  a 
critical  hydrogeologic  threshold;  thus,  the  Yukon  Flats  Basin  may  be  on  the  verge  of  major 
hydrologic  change  should  the  current  permafrost  distribution  degrade.  This  possibility 
underscores  the  need  for  improved  characterization  of  permafrost  and  other  hydrogeologic 
information  in  the  study  area  via  geophysical  techniques  and  ground-based  observations.  The 
regional  groundwater  flow  model  identifies  the  potential  importance  of  lake  taliks  in  enhancing 
groundwater  circulation  in  the  system. 

Evolution  of  and  flow  through  lake  taliks  has  been  explored  in  detail  using  SUTRA  with  new 
freeze/thaw  capabilities  that  capture  both  rates  of  thaw  and  thaw  patterns.  Results  confirm  the 
importance  of  simulating  advective  heat  transport  in  addition  to  conduction  as  an  important 
component  in  accelerating  thaw  rates.  The  response  of  lake  systems  to  lake  area  decline  has  been 
explored  using  both  analytical  and  SUTRA  simulation  approaches.  The  shallow,  continuous 
gravel  layer  in  the  Yukon  Flats  is  particularly  vulnerable  to  permafrost  degradation,  potentially 
resulting  in  large  changes  in  watershed  runoff,  shallow  groundwater  flow,  and  hydrologic 
connectivity  between  lakes  and  stream  networks.  Modeling  results,  supported  by  the  geophysical 
field  investigations,  indicate  that  permafrost  aggradation  in  dried  lake  margins  may  offer  some 
resilience  to  lake  area  decline,  but  favorable  conditions  for  the  reformation  of  permafrost  may  be 
countered  by  climate  warming.  Cross-scale  investigations  integrating  characterization  and 
dynamic  cryo-hydrologic  processes  are  needed  to  further  constrain  the  fate  of  lake  and  wetland 
distribution  in  the  Yukon  Flats  and  other  arctic/subarctic  lowlands. 


6,2  Implementation  of  Tools  and  Approaches  Developed 

6.2.1  SUTRA 

The  uses  SUTRA  code,  as  a  result  of  development  and  enhancement  work  done  under  this 
project,  has  new  general  capabilities  that  will  allow  scientists  and  practitioners  worldwide  to 
evaluate  ground-ice  and  permafrost  relations  within  hydrologic  systems,  and  thereby  to  evaluate 
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impacts  of  climate  changes  on  earth-surfaee  environments.  As  with  all  USGS  publie  hydrologie 
software,  the  new  versions  of  the  SUTRA  eode  (both  programs  and  exeeutables)  are  made 
available  online,  free  of  eharge,  to  everyone.  In  addition,  the  eode  doeumentations  (i.e.  user 
manuals)  and  simulation  preproeessing  software  are  made  available  on  the  same  website  for  free 
download.  Thus,  all  USGS-produeed  eomputer  eode  and  information,  in  the  present  ease  for 
SUTRA,  is  published  online  in  open  aecess  format.  The  US  publie,  government  ageneies,  DoD 
interests,  universities  and  private  interests,  as  well  as  similar  interests  in  other  eountries, 
regularly  aeeess  the  website  for  hydrologie  information  and  downloading  software. 

As  with  previous  new  versions  of  the  USGS-SUTRA  eode,  the  enhaneed  versions  produeed 
under  this  SERDP  projeet  will  eontinue  to  be  taught  in  USGS-led  eourses  and  workshops. 
Teehnology  transfer  is  part  and  pareel  of  all  new  versions  of  SUTRA  released  by  USGS.  As  part 
of  this  projeet,  two  workshops  have  already  been  taught  that  foeus  on  the  new  freeze/thaw 
eapabilities  of  SUTRA  (at  USGS  in  Lakewood  CO  and  at  MeGill  University,  loeation  of 
Professor  J.  MeKenzie,  one  of  the  eode  developers  and  SERDP  key  performers,  in  Montreal, 
Canada).  A  number  of  presentations  about  the  eode  and  about  results  it  has  already  given  have 
been  made  during  the  duration  of  this  SERDP  projeet  (domestie  and  foreign  universities  and 
institutes,  and  at  USA  national  meetings  such  as  AGU  and  GSA),  and  more  are  planned  by  all 
investigators.  In  2015,  Principle  Investigator,  C.  Voss,  is  providing  a  leeture  about  cold-regions 
freeze/thaw  hydrology  with  a  foeus  on  results  of  this  projeet,  partieularly  those  obtained  with  the 
new  SUTRA  eode,  at  universities  and  institutes  throughout  the  USA  and  in  some  other  eountries. 
This  is  part  of  the  Birdsall-Dreiss  Leetureship,  sponsored  by  the  Geologieal  Soeiety  of  America 
(http://gsahvdrogeologv.org/BirdsallDreiss2015.htm) ,  for  whieh  Voss  has  been  designated  as  the 
2015  lecturer.  Voss  expeets  to  present  30  to  50  leetures  during  2015,  a  significant  technology 
transfer  effort  for  the  hydrologie  community  that  will  make  many  professionals  and  students 
aware  of  progress  made  under  the  SERDP  projeet. 

Moreover,  key  faetors  in  transitioning  SUTRA  to  the  user  eommunity  have  always  been  (1) 
USGS  partieipation  in  and  support  of  new  projeets  initiated  by  users  of  SUTRA,  and,  (2)  reliable 
user  support  based  on  personal  eontact  via  telephone  and  email  of  SUTRA  users  worldwide.  In 
the  ease  of  SUTRA  freeze/thaw  eapabilities,  user  support  regarding  eode  use  and  problems 
eneountered  is  being  provided  by  key  performers  Professor  J.  MeKenzie  (MeGill  University)  and 
A.  Provost  and  by  PI  C.  Voss  (USGS).  Eor  questions  related  to  the  simulation  preproeessing  and 
post-proeessing  software,  R.  Winston  (USGS)  provides  support.  As  part  of  eode  development,  a 
beta  version  of  SUTRA  with  freeze/thaw  eapabilities  had  been  distributed  to  seleeted  beta  testers 
during  the  SERDP  projeet  period,  and  together  with  the  SERDP-funded  eode  developers  at 
USGS,  these  beta  testers  have  applied  the  code  to  study  a  wide  variety  of  eold-regions  settings. 
Results  of  these  studies  have  already  been  published.  Thus,  transitioning  the  eode  to  the 
eommunity  is  a  eontinuous  proeess  that  oeeurred  simultaneously  with  eode  development  and 
enhaneements.  A  small  but  signifieant  user  base  has  thus  already  been  established  for  SUTRA 
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with  freeze/thaw  capabilities,  with  most  users  outside  of  USGS.  This  user  base  will  eontinue  to 
grow. 

6.2.2  Geophysieal  Methods  for  Cold  Regions 

Our  multi-seale  geophysieal  approaeh  is  being  transferred  through  published  and  planned  peer- 
review  publieations,  as  well  as  direct  training  to  cooperators  in  Alaska.  In  August  2014,  project 
Pi’s  delivered  a  3-day  training  course  at  University  of  Alaska,  Fairbanks,  with  participants  from 
academia,  federal  and  state  ageneies,  and  the  private  seetor.  Course  materials  will  be  made 
available  online  following  USGS  review  and  approval. 

A  geophysical  software  package  for  frequency-domain  electromagnetic  inversion  (FEMIC)  was 
developed  and  demonstrated  under  this  project.  The  development  of  this  software  fills  an 
important  gap  in  geophysical  analysis  since  eommereially  available  software  is  not  amenable  to 
inversion  of  large,  aerially  distributed  eleetromagnetie  data,  sueh  as  eolleeted  for  this  projeet. 
The  FEMIC  software  is  eurrently  undergoing  internal  USGS  testing  prior  to  publieation.  The 
eode  and  exeeutable  will  be  publieally  available,  and  will  have  a  broad  range  of  applieation 
ineluding  mapping  of  shallow  permafrost  via  frequeney  domain  eleetromagnetie  surveys. 

Understanding  the  hydrogeophysieal  responses  to  permafrost  dynamies  under  different 
hydrologie  and  elimatie  eonditions,  and  in  different  geologieal  settings,  is  important  for  guiding 
the  interpretation  of  existing  geophysieal  datasets  and  also  for  planning  future  surveys. 
Geophysieal  models  are  inherently  uneertain  and  ambiguous  beeause  of  (1)  the  resolution 
limitations  of  any  geophysieal  method  and  (2)  the  weak  or  non-unique  relationship  between 
hydrologie  properties  and  geophysieal  properties.  The  general  framework  for  eoupling 
geophysieal  predietions  to  hydrologic  simulations  of  permafrost  evolution,  ineluding  a  novel 
physieal  property  relationship  that  aeeounts  for  the  eleetrical  response  to  changes  in  lithology, 
temperature,  and  iee  eontent,  as  well  as  a  rigorous  analysis  of  geophysieal  uncertainty,  is  readily 
transferrable  to  other  eold-region  studies  of  geophysieal  charaeterization.  Although  the  foeus 
here  was  on  AEM  data,  other  types  of  eleetrieal  or  eleetromagnetie  measurements  eould  be 
readily  simulated  using  the  same  resistivity  model.  Eor  example,  future  efforts  could  focus  on 
the  simulation  of  other  types  of  geophysical  data  (e.g.  nuclear  magnetie  resonanee  or  ground 
penetrating  radar)  using  the  same  basie  modeling  approaeh. 

6.2.3  Integrated  Approaeh  for  Euture  Assessments 

A  eomprehensive  charaeterization  of  lateral  and  vertieal  distribution  of  permafrost  may  require 
complementary  ground-based  and  aerial  geophysics,  developed  as  part  of  this  projeet,  depending 
on  the  seale  of  interest.  Knowledge  of  permafrost  distribution  is  a  fundamental  underpinning  of 
effeetive  management  and  possible  mitigation  of  elimate-ehange  impaets  of  ground-iee  changes 
on  hydrology  and  land-surfaee  environments.  It  is  neeessary  to  know  where  there  is  ground  iee  in 
order  to  prediet  the  ehanges  in  hydrology  and  land  surfaee  environment  that  may  oeeur  should 
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the  ground  ice  disappear  or  modify  its  spatial  distribution.  Aerial  geophysics  provides  ground-ice 
distributions  over  distances  of  tens  of  kilometers,  whereas  ground-based  geophysics  provides 
more  details  of  ground-ice  distribution  over  distances  of  hundreds  of  meters.  Linkage  and 
correlation  of  geophysical  characterizations  with  satellite  remote  sensing  data  may  allow  for 
generalization  of  shallow  ground-ice  distribution  to  regions  not  yet  evaluated  by  geophysics,  but 
such  upscaling  techniques  require  further  development  and  additional  ground  truth  information. 

The  new  numerical  simulation  tool,  SUTRA  with  groundwater  freezing  and  thawing,  allows 
scientific  hypothesis  testing  to  be  carried  out  to  define  the  most  important  factors  that  control 
subsurface  hydrologic  and  thermal  response  to  conditions  on  the  ground  surface  in  a  particular 
setting.  Ground  surface  conditions  that  can  be  considered  include,  insolation,  vegetation  and 
shading,  vegetation  loss  or  regrowth,  air  temperature,  surface  slope,  surface  hydrology  (bare  or 
surface  water),  organic  vs  mineral  soil  layers,  fires,  and  human  infrastructure  (e.g.  roads, 
drainage  channels,  structures).  Elucidation  of  rates  of  change  and  key  factors  controlling  the 
processes  of  concern  in  managing  various  settings  in  DoD  lands,  for  example  using  SUTRA, 
must  be  accomplished  before  effective  management  strategies  can  be  defined. 

An  example  that  combines  the  above  approaches  concerns  the  fate  of  a  lowland  lake  in  Alaska, 
studied  as  part  of  this  project.  A  basic  hydrologic  question  concerns  the  evolution  of  Twelvemile 
Lake’s  water  level,  which  had  been  dropping  steadily  for  30  years.  Disappearance  of  lakes  is  a 
widespread  occurrence  in  interior  Alaska.  Airborne  and  ground-based  geophysics  allowed 
determination  of  the  ground-ice  and  permafrost  distribution,  including  newly  developed 
permafrost,  around  the  lake.  Hydrologic  evaluation  of  the  lake’s  water  budget,  using  basic 
quantitative  analysis  and  using  simulation  with  SUTRA,  allowed  identification  of  the  candidate 
factors  that  can  control  lake  water  level.  Further  analysis,  conducted  after  sudden  flooding  and 
refilling  of  the  lake  in  2013,  found  that  the  location  of  the  permafrost  table  can  be  an  important 
control  on  hydraulic  connectivity  between  upstream  lakes  and  Twelvemile  Lake.  Depending  on 
the  permafrost  distribution,  upstream  lakes  may  be  able  to  refill  Twelvemile  Lake  during  Yukon 
River  flood  events,  and  may  subsequently  control  the  amount  of  Twelvemile  Lake  water  gain 
and  loss  on  a  yearly  basis.  This,  in  combination  with  the  amount  of  snow  that  falls  in  the 
Twelvemile  Lake  watershed  controls  the  rate  of  decrease  in  lake  level  in  the  years  after  a 
refdling  event.  Snowfall  depends  on  climate  and  snow  capture  on  the  ground  surface  depends  on 
the  type  of  vegetation  that  exists  in  the  watershed.  Thus,  this  particular  system  is  an  example  of 
how  the  functioning  of  a  particular  setting  can  be  elucidated  using  the  methods  developed  in  this 
project.  As  an  example  of  an  intervention,  consider  the  possibility  that  this  watershed  is  a  DoD 
concern  in  which  it  was  desired  that  the  lake  level  remain  high,  rather  than  drop  during  the  years 
between  flood  events.  Perhaps  a  re-vegetation  program  could  be  implemented  after  geophysical 
characterization  that  indicates  the  location  of  permafrost  and  after  evaluation  of  potential 
intervention  approaches  using  SUTRA  simulation.  The  program  might  then  consist  of:  (1) 
particular  regional  vegetation  would  be  planted  that  captures  snow,  providing  a  greater  source  of 
spring  meltwater  to  the  lake,  and,  (2)  the  vegetation  would  shade  the  ground  surface,  keeping  the 
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permafrost  table  shallow  (this  would  quickly  route  meltwater  to  the  lake,  with  less 
evapotranspiration  losses  that  would  occur  should  the  permafrost  table  be  deeper). 

6,3  Ramifications  for  DoD  Operations,  Infrastructure,  and  Planning  in  Interior  Alaska 

The  results  and  products  derived  from  this  study  can  help  support  decision  making  for  DoD 
infrastructure  and  training  activities  in  interior  Alaska,  particularly  in  areas  most  vulnerable  to 
climate  change  and  land  use  disturbance  via  permafrost  degradation  and  resultant  changes  in 
hydrogeology.  The  U.S.  Army  is  the  largest  DoD  land  user  in  Alaska  and  conducts  all  season 
training  activities  throughout  its  1.5  million  acre  training  areas  in  the  Alaskan  Interior  (Figure 
6.3.1).  Most  of  the  impact  ranges  used  by  the  Air  Force  in  Alaska  are  located  on  Army  lands,  and 
the  Army  and  Air  Force  stage  joint  training  exercises  on  Army  lands.  The  main  DoD  sites  are 
located  between  Fairbanks  and  Delta  Junction  and  include  Fort  Wainwright  and  its  associated 
training  areas  (Tanana  Flats,  Donnelly  East  and  Donnelly  West),  Eielson  Air  Force  Base,  and 
Fort  Greely. 


Figure  6,3,1  An  elevation  map  of  interior  Alaska  DoD  training  areas  and  cantonments. 
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The  cantonments  of  Fort  Wainwright  and  Eielson  Air  Force  Base,  the  Tanana  Flats  Training 
Area,  and  most  of  Donnelly  Training  Area  West  and  East  are  located  on  lowland  landscapes 
underlain  by  discontinuous  ice-rich  permafrost.  This  complex  mosaic  of  permafrost  and 
permafrost-free  land  is  predominantly  low  gradient  (1-2  m  km"')  and  even  small  changes  in 
surface  topography,  vegetation,  and  the  permafrost  thermal  regime  can  have  dramatic  effects  on 
surface  and  shallow  subsurface  hydrology. 

Quantitative  predictive  models  are  needed  for  understanding  the  processes  occurring  in  the 
shallow  subsurface  in  cold  regions  in  response  to  infrastructure,  for  predicting  possible  future 
changes  exacerbated  by  climate  change,  and  for  aiding  decision  making  related  to  infrastructure. 
Approaches  and  tools  developed  in  this  study  could  be  applied  to  help  inform  DoD  decision 
making,  particularly  where  numerical  modeling  combined  with  geophysical  characterization 
could  be  used  to  evaluate  the  vulnerability  (thaw  settlement,  ponding  etc..)  of  areas  to 
groundwater/permafrost  interactions  and  project  likely  consequences.  Such  information  is 
essential  for  accurately  assessing  the  structural  integrity  of  both  built  and  proposed 
infrastructure.  Eor  example,  the  evolution  of  permafrost  degradation  affected  by  and  affecting 
groundwater  flow  could  be  simulated  to  help  guide  the  siting  or  design  of  roads,  runways  and 
airfields,  buildings,  bridges,  and  supportive  engineering  structures.  Talik  development  and 
enhanced  subsurface  hydrologic  connectivity  in  response  to  climate  change  and  other 
disturbances  including  infrastructure  development  could  greatly  alter  local  hydrology, 
particularly  near  wetlands.  A  major  challenge  for  road  construction  is  in  understanding  how  and 
where  earthen  berms  will  alter  local  and  regional  surface  water  hydrology  in  the  low  gradient 
terrain.  Physics-based  modeling,  particularly  using  SEITRA  and  the  enhancements  made  within 
this  project,  could  help  elucidate  permafrost-controlled  surface-water/groundwater  interactions. 
The  fen  and  bog  features  in  lowland  terrains  are  dynamic  wetland  landscapes  that  are  also 
intimately  related  to  the  presence  (or  absence)  of  permafrost.  They  provide  natural  fire  breaks 
and  support  numerous  species  of  interest  for  natural  resources  planning.  However,  with  the 
projected  climate  warming  in  the  area  (5°C  by  2100;  Chapman  and  Walsh,  2007)  many  of  these 
features  could  be  altered. 

The  results  from  this  effort  and  from  the  team  representing  SERDP  Project  RC-2110 
(Addressing  the  impacts  of  Climate  Change  on  US  Army  Alaska  with  Decision  Support  Tools 
Developed  through  Eield  Work  and  Modeling)  could  be  combined  to  provide  planning  level 
information  for  DoD  in  interior  Alaska.  The  measurements  from  fire  disturbed  landscapes  and 
fens  in  the  Tanana  Plats  lowlands  could  provide  critical  input  to  a  suite  of  SUTRA  models 
designed  to  examine  current  and  projected  hydrologic  surface/subsurface  connectivity  and 
implications  for  surface  water  distribution  in  that  region.  Eor  upland  regions,  like  the  Yukon 
Training  Area  and  parts  of  Donnelly  Training  Areas  East  and  West,  SUTRA  application  to 
representative  hillslopes  could  be  used  to  identify  the  thermal  and  hydrologic  response  (and 
evaluate  resilience)  to  forest  fires  and  to  provide  critical  information  on  current  and  projected 
summer  season  flows  in  watersheds  in  the  area. 
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Information  on  the  sizing  of  culverts  and  bridges  and  the  siting  of  winter  stream  crossings  could 
be  gleaned  from  modeling  efforts  designed  to  address  projected  increased  stream  baseflow  across 
discontinuous  permafrost  landscapes  with  continued  climate  warming.  This  is  of  particular 
interest  for  anadromous  streams  where  a  higher  level  of  scrutiny  and  more  robust  engineering 
designs  are  required  to  ensure  fish  vitality  and  health.  Many  of  the  streams  in  Xanana  Flats  and 
Donnelly  West  have  never  been  gaged  and  they  have  strong  seasonal  variations  in  flow.  Some 
streams  have  glacial  sources  and  others  are  predominantly  spring  fed.  However,  little  is 
understood  about  the  role  of  permafrost  on  stream  flows  in  these  systems.  In  locations  where 
permafrost  degradation  contributed  to  enhanced  groundwater  flow  and  stream  baseflow,  the 
changes  in  the  seasonal  hydrograph  and  stream  temperatures  could  be  dramatic.  Both  factors 
strongly  impact  fish  habitat. 

The  ground  based  geophysical  methods  used  to  map  permafrost  distribution,  characterize  annual 
ground  freeze/thaw  cycles,  and  monitor  ground  ice  changes  at  our  sites  could  be  applied  to  a 
variety  of  settings  to  support  the  siting  of  buildings  and  could  be  combined  with  modeling  efforts 
to  help  predict  the  likely  response  of  permafrost  to  climate  and  anthropogenic  disturbance.  When 
coupled  with  airborne  geophysical  measurements  the  ground  based  measurement  results  could  be 
more  broadly  applied  to  the  remote  locations  in  training  ranges  where  no  road  access  is  currently 
available.  There  are  instances  where  infrastructure  in  poorly  characterized  areas  has  been  over¬ 
engineered,  designed  for  thaw-unstable  permafrost  in  an  area  where  permafrost  was  thaw-stable. 
A  carefully  designed  approach  integrating  geophysical  and  permafrost/hydrology  techniques 
developed  in  our  study  could  help  minimize  unnecessary  costs  and  wisely  allocate  resources 
associated  with  proposed  infrastructure  projects. 


6,4  Remaining  Gaps  in  Knowledge 

Geophysical  characterization  performed  under  the  auspices  of  this  project  was  limited  in  spatial 
extent.  Large  gaps  in  the  spatial  extent  and  depth  of  permafrost  still  exist  over  many  areas  of 
interior  Alaska.  The  airborne  and  ground-based  geophysical  methods  explored  in  this  project 
offer  a  comprehensive  overview  of  current  techniques.  Analysis  of  these  methods  offers 
guidance  on  the  tradeoffs  between  resolution,  depth  of  interrogation,  and  spatial  coverage.  Such 
guidance  will  be  useful  for  designing  future  permafrost  characterization  campaigns  to  fill 
important  baseline  gaps  in  interior  Alaska. 

Recent  developments  in  numerical  models  to  simulate  coupled  heat  and  water  transport  in 
temperate  to  freezing  subsurface  conditions  have  advanced  our  ability  to  understand  water 
movement  in  permafrost  environments.  However,  several  gaps  exist  in  understanding  and 
physically  describing  basic  processes  and  property  relationships  associated  with  freeze/thaw 
cycles,  particularly  for  variably  saturated  conditions.  Accepted  practices  for  implementing  the 
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mechanisms  of  cryosuction,  frost  heaving,  and  thermokarsting  (land  subsidence  associated  with 
shallow  permafrost  thaw  and  soil  matrix  collapse)  in  physics-based  models  have  not  been  widely 
established;  these  are  areas  of  active  research  and  development.  Also,  the  responses  of  hydraulic, 
thermal,  and  mechanical  properties  of  soil  across  the  freeze/thaw  transition  for  various  soil  types 
are  not  well  known.  Very  limited  attention  has  been  given  to  soil  freezing  properties  that 
influence  the  flow/impedance  of  water  and  mechanical  strength  of  soil  under  variably  saturated 
conditions  in  the  subsurface.  Additional  laboratory  and  field  investigations  are  needed  to  better 
constrain  thermal,  mechanical,  and  hydraulic  relationships  for  representation  in  cold-region 
model  simulations. 

Understanding  connections  between  thermally-influenced  landscape  processes  and  inland  waters 
(rivers,  streams,  lakes,  wetlands)  in  permafrost  systems  represents  a  broad  area  where  additional 
research  is  needed.  For  example,  hydrologic  responses  to  land  subsidence,  erosion,  and  changes 
in  channel  geomorphology  derived  from  permafrost  thaw  and  climate  warming  are  not  well 
quantified  nor  are  potential  thermal  feedbacks  to  these  landscape  processes.  Improved 
knowledge  of  groundwater-surface  water  interaction  in  permafrost  systems  is  required  for 
predicting  the  trajectories  of  inland  water  bodies  under  a  warming  climate.  Approaches  for 
addressing  these  critical  gaps  will  require  a  combination  of  remote  sensing,  geophysical  and 
hydrologic  field  characterization,  and  coupled  thermal-hydrologic  modeling. 


6,5  Recommendations  for  Long-Term  Permafrost  Monitoring  Program 

Multi-method,  multi-scale  approaches  to  characterize  permafrost  are  needed  to  capture  dynamic 
processes.  Application  of  many  geophysical  methods  used  and  developed  in  RC-2111  (and  also 
distributed  temperature  sensor  cable  employment)  may  be  used  to  monitor  thermal  and  physical 
state  of  permafrost  at  extensive  sites  in  different  landscapes  and  for  calibration  data  for 
airborne/satellite  platforms.  Several  of  the  geophysical  instruments  used  here,  along  with  thermal 
imaging,  are  amenable  to  future  deployment  as  payloads  on  unmanned  aerial  vehicles 
(unmanned  aerial  vehicles,  i.e.,  drones),  enabling  cost-effective  mapping  in  remote  areas.  Long¬ 
term  shallow  temperature  monitoring  would  be  useful  for  defining  appropriate  boundary 
conditions  in  SUTRA  and  other  model  applications.  Extensive  permafrost  characterization  and 
monitoring  sites  should  span  a  spectrum  of  vegetation  types  and  fire  history  in  order  to  provide  a 
comprehensive  assessment.  Sites  could  be  collocated  with  military  lands  or  areas  of  military 
interest  to  enhance  DoD  relevance.  Monitoring  equipment  installed  at  Beaver  Meadow  Lake  and 
Twelvemile  Lake  as  part  of  the  RC-2111  project  continues  to  collect  pressure  and  temperature 
data.  Luture  projects  could  benefit  from  this  ongoing  data  collection. 

Monitoring  of  physical  changes  in  the  landscape  that  are  reflective  of  permafrost  thaw,  such  as 
land  subsidence,  thermokarst,  and  erosion  should  be  included  in  a  permafrost  monitoring  plan. 
Airborne  approaches  (i.e.,  LiDAR,  interpretation  of  thermokarst  feature  using  high-resolution 
satellite  imagery  and/or  aerial  photos,  etc...)  may  be  well  suited  for  monitoring  geomorphic 
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change  over  large  seales.  In  addition,  basie  surfaee  and  subsurfaee  hydrology  data  should  be 
ineluded  in  a  monitoring  plan  to  understand  implieations  assoeiated  with  permafrost  ehange. 
Sueh  data  inelude,  but  are  not  limited  to  soil  moisture,  water  table  elevation,  lateral  and  vertieal 
hydraube  head  distribution,  rain  and  snowfall,  evaporation,  runoff,  seasonal  streamflow,  and  lake 
level.  Basie  surfaee  and  subsurfaee  hydrology  data  should  be  ineluded  in  permafrost  monitoring 
to  understand  implieations  assoeiated  with  permafrost  ehange.  Sueh  data  inelude,  but  are  not 
limited  to  soil  moisture,  water  table  elevation,  lateral  and  vertieal  hydraube  head  distribution, 
rain  and  snowfall,  evaporation,  runoff,  seasonal  streamflow,  and  lake  level. 

Means  to  ground-truth  and  verify  interpretations  are  fundamental  to  any  long-term  monitoring 
plan.  We  suggest  inelusion  of  sentinel  deep  wells  (for  hydrologie,  thermal,  geoehemieal,  and 
geophysieal  monitoring)  loeated  in  areas  that  will  provide  most  information  on  permafrost  and 
hydrologie  evolution,  in  sensitive  eeosystems  and  in  areas  with  threatened  infrastrueture. 

Modeling  studies  using  the  approaehes  developed  in  RC-2111  ean  be  used  to  guide  a  strategie 
permafrost  monitoring  plan  in  (1)  identifying  areas  of  high  vulnerability,  (2)  eonstraining 
appropriate  measurement  frequeney,  and  (3)  interpreting  monitoring  results. 

Coordination  with  other  existing  and  planned  monitoring  programs  (e.g.  Department  of  Energy’s 
Next-Generation  Eeosystem  Experiments,  NASA’s  Aretie  Boreal  Vulnerability  Experiment, 
National  Seienee  Eoundation’s  Aretie  Observing  Network,  USGS’s  Alaska  Climate  Seienee 
Center,  and  National  Snow  &  lee  Data  Center’s  Cireumpolar  Aetive  Eayer  Monitoring  Network) 
will  alleviate  required  resourees  and  maximize  information. 
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